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During this period investigations were carried out in two areas: 

1. Energy- and structure-related properties of small gold clusters 
deposited on the GaAs(llO) surface were investigated in this 
work using a molecular dynamics procedure. A recently devel- 
oped potential energy function based on two- and three-body 
interactions was employed in calculating energies and forces. 
These calculations produced some consistent results with ex- 
periments. Calculations indicate that Au atoms adsorbed on 
the GaAs(llO) surface do not favor the formation of 3D clus- 
ters at the very early stages of deposition (i.e, up to Au 5 ). In all 
cases, gold atoms were found to prefer sites near Ga atoms. Sites 
for Au near As atoms are energetically less favorable. Further- 
more, the present study suggests that three-body interactions 
involving triplets of Au and As atoms play an important role in 
determining sites and binding energies for deposited Au atoms 
on the GaAs(llO) surface. Results and the method of calcula- 
tion are given in Appendix 1. 

2. A comparative study was conducted in this part for six clas- 
sical many-body potentials developed recently for silicon sys- 
tems. Extensive static calculations were performed using these 
potentials on small Si clusters, bulk point defects, elastic con- 
stants, polytypes, pressure induced phase transformations and 
low index plane surfaces. Similarities and differences between 
six potentials were identified and their transferability as well as 
their accuracy with respect to experiment and first-principles 
methods were assessed. In general, all these potentials do a 
relatively poor job of modeling the energetic of small clusters 
as well as the various reconstructions of the Si(lll) surface. 
They, however, provide a fair to good description of the prop- 
erties of bulk diamond cubic silicon, its intrinsic defects, and of 
the Si(100) surface. Besides the fact that none of them model 
7 r-bonding, their inability to be more transferable lies in their 
inadequate description of the angular forces. Each potential has 
its strengths and limitations but none appears to be clearly su- 
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perior to the others and none is totally transf err able, However, 
despite their shortcomings we feel that some of these potentials 
will be useful in large scale simulations of materials-related prob- 
lems. They can give valuable insights into phenomana which are 
otherwise intractable to investigate either experimentally or via 
first-principles methods. Details of this investigation are pre- 
sented in Appendix 2. 
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SIMULATION CALCULATIONS FOR GOLD CLUSTERS 
ON THE GaAs(llO) SURFACE 
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Introduction 


Due to its increasing technological importance in recent years, 
numerous experimental as well as theoretical investigations have been 
carried out in the area of metal-semicontactor interfaces. Today, 
structural aspects of such interfaces related to atomic configurations 
in particular, are highly desired. Simulation calculations are carried 
out here to analyze some of the energy- and structure-related proper- 
ties of the Au/GaAs interface at an atomic level. In this study, small 
clusters of gold with varying number of atoms were deposited on the 
GaAs(llO) surface. Calculations were carried out to determine bind- 
ing energies of the clusters as well as the energetically most favorable 
binding sites for deposited gold atoms. 


Method of Calculation 


A molecular dynamics procedure based on the Nordsieck-Gear al- 
gorithm was employed here. Calculations were performed considering 
a potential function based on two- and three-body interactions. The 
total energy of interaction of a system of N particles was calculated 

aS ‘ jv N 

$ = u{fj , Tj) + 'y ^ u(fi,Tj,r k ,) ( 1 ) 

i<j i<j<k 

where, and «(?<,*>-, f»), denote the two- and three-body in- 

teractions, respectively. The two-body part was represented by the 
Lennard- Jones pair potential: 
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where, r v = |f; - fj |; r 0 represents the equilibrium distance and e 
denotes the two-body energy at ry - r 0 . For the three-body part the 
Axilrod-Teller triple dipole potential was taken into consideration: 
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where, »i,Sj,»k and r ih r ik ,r jk represent the angles and the sides of 
the triangle formed by the three particles i, j and k. respectively. The 
three-body intensity parameter is denoted by Z . 

It has been demonstrated that this potential function is able to 
reproduce several energy- and structure-related properties for GaAs 
systems [1-4]. In a recent study this same potential function has been 
parametrized for the Au-Ga-As system and, at the same time, it has 
been used for determining high energy adsorption sites for a single 
gold atom deposited on the (110) surface of GaAs [5]. Parameters of 
the potential function for the Au-Ga-As system are given in Table 1. 
In the present investigation a similar procedure was employed and cal- 
culations were carried out for sm a ll clusters of gold atoms (containing 
up to five Au atoms) deposited on the (110) surface of GaAs. Simu- 
lations were performed for the low temperature case (T«l K), and a 
time step of 1.4 x 10 -15 sec. was considered throughout this study. 
Also, periodic boundary conditions were employed in two directions 
to provide continuity for the exposed surface. The (110) surface was 
generated, first, as an abrupt termination of a properly oriented GaAs 
lattice. Next, Au atoms were deposited on this ideal surface and the 
system was equilibrated under the molecular dynamics code. In gen- 
eral, for varying numbers and positions of gold adatoms deposited on 
GaAs(llO), an average of 8,000 molecular dynamics steps was found 
to be sufficient for a complete equilibration of the system. In the 
present work, however, to assure that the system is fully equilibrated, 
10,000 time steps were employed. For the very last 2,000 time steps 
the temperature rescaling was turned off and calculated values, in 
each case, were collected for averaging. During this period the tem- 
perature remained fairly constant indicating that the system was fully 

relaxed. 

Calculations for the (110) surface were carried out employing a 
GaAs substrate of 60 unit cells in a (4x3x5) arrangement containing a 
total of 240 atoms. The exposed surface consisted of 12 (4x3) surface 
cells. A portion of this surface is depicted in Figure 1. In addition to 
Ga and As atom positions in the top two layers of the substrate, vari- 
ous energetically favorable sites for Au adatoms are also shown. In the 
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calculations, deposited Au atoms as well as substrate atoms located 
in the top two layers were permitted to fully relax. The rest of the 
substrate atoms in the system were held fixed in their lattice sites, 
but they were permitted to fully interact with the relaxing atoms. 
Because the calculations in this study were carried out at a low tem- 
perature, Au atoms positioned at the surface often remain trapped at 
the nearest minimum. To eliminate this undesired situation, there- 
fore, calculations for each cluster, were repeated considering several 
different initial configurations and only those with lowest energies are 
reported here. 


Results and Discussions 

Figure 1 shows favorable binding sites for Au atoms of different 
clusters deposited on the (110) surface of GaAs. For a single Au atom, 
one of the favorable binding sites is shown by the letter ‘a’ m Figure 1. 
The location of this site along with the binding energy value given m 
Table 2 for a single Au atom, are consistent with our earlier report [5]. 
In the case of Au 2 , the most favorable sites for Au atoms are ‘a’ and 
‘b’ positions which are identical high energy sites in two neighboring 
surface cells. In this case the separation between two Au atoms is 
quite large (~ 4 A). Therefore, for a dimer deposited on the (110) 
surface the Au-Au interaction is expected to be rather weak and this 
is reflected in the e b value given in Table 2. For Au 3 , sites ‘a’, ‘b 
and ‘c* were predicted to be most favorable. The deposited cluster, 
in this case, is an isoceles triangle and the distance between two Au 
atoms (forming the equal sides of the triangle) is about 2.5 A. This 
is slightly shorter than the Au-Au equilibrium distance (see Table 1). 
Even though pair-wise interactions between gold atoms are expected 
to be more attractive, the Au atom at site V is situated in close 
proximity to three As atoms (two are in the top layer and the other 
is in the second layer) which form Au-Au- As and Au-As-As trimers 
with acute angles providing strong repulsive three-body forces. This is 
reflected in the e^ value for Au 3 which is about 0.21 eV/per Au atom 
weaker in binding energy than the individual gold atom (see Table 2). 
For Au 4 , sites ‘a’, ‘c’, c d’ and *e’ were found to be the most favorable 
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adsorption sites. In this case, while three Au atoms are located near 
each other, the forth one, which is in site ‘e\ is somewhat separated. It 
is equidistant from atoms at ‘c’ and ( d’ with a separation of about 4.7 
A. At that separation no significant Au-Au interactions are expected 
coming from the Au atom positioned at the site V. Repulsive three- 
body interactions, however, in this configuration, arising from Au- 
Au-As and Au-As-As triplets, were reduced. Accordingly, the value 
of ej, for Au 4 was found to be decreased indicating that the bonding 
is somewhat stronger in this case. For Au 5 , favorable adsorption 
sites were found to be ‘a’, ‘c\ ‘d\ T and ‘g’ as indicated in Figure 

1. In this case the two Au atoms located at sites T and ‘g’ are 
somewhat separated from atoms at ‘a’, ‘c’ and c d\ Like Au 4 , in this 
case also farther atoms (located at sites T and ‘g’) are not expected to 
contribute much Au-Au interaction. Similarly, calculations indicate 
that repulsive three-body interactions (due to triplets involving Au 
and As atoms) also remained at a low level. Therefore, the value of 
the binding energy per atom, e&, for the Au 5 case as shown in Table 

2, did not change much when compared with Au 4 . 

Based on the potential energy parameters given in Table 1, it 
is expected that an isolated Au-Ga dimer would be less energetic 
than an isolated Au-As dimer. On the GaAs (110) surface, how- 
ever, Au adatoms in general, were found to prefer sites closer to Ga 
atoms. Present calculations clearly indicate that this is because of 
strong three-body interactions exhibited by trimers involving Au and 
As atoms. Judged by the values of the three-body parameters given 
in Table 1, on the other hand, energies coming from trimers involving 
Au and Ga atoms are expected to be rather small. The three-body 
part of the potential function, as given in Eq. 3, provides strong repul- 
sive energies, in particular, for Au-Au-As and Au-As-As trimers with 
acute angles. Accordingly, surface sites for Au near As atoms become 
energetically less favorable. This outcome is in general qualitative 
agreement with various reports [6,7] indicating that Au is prefentially 
bonded to Ga atoms at the surface. Feenstra [6] using an STM tech- 
nique measured the lateral distance between a Au adatom and a top 
layer substrate Ga atom as ~1.4 A. While this value is consistent 
with the calculated lateral distance of 1.44 A between the single Au 


5 



adatom and its closest Ga neighbor, present calculations predict that 
this nearest neighbor Ga atom is located in the second layer. 


Conclusions 

Calculations indicate that Au atoms adsorbed on the GaAs(llO) 
surface do not favor the formation of 3D clusters at the very early 
stages of deposition (i.e, up to Au 5 ). In all cases, gold atoms were 
found to prefer sites near Ga atoms. Furthermore, the present study 
suggests that three-body interactions involving triplets of Au and As 
atoms play an important role in determining sites and binding en- 
ergies for deposited Au atoms on the GaAs(llO) surface. Results 
obtained in this investigation are strictly based on energetics and no 
entropic aspects were taken into consideration here. Therefore, it is 
recommended that extreme care should be exercised when comparing 
these results with experimental findings. 
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Table 1. Parameters for the Potential Energy Function [2]. 




Two-Body Part 


€ (eV) 

?o (A) 

Au — Au 

0.9760 

2.6685 

Ga — Gcl 

1.0039 

2.4607 

As — As 

1.1641 

2.4913 

Au — Ga 

0.8860 

2.5540 

Au — As 

1.7500 

2.5350 

Ga — As 

1.7379 

2.4481 



Three-Body Part 



Z (eV A 9 ) 

Au - Au - Au 


2009.0 

Ga-Ga- Ga 


1826.4 

As - As - As 


2151.9 

Au — Au — Ga 


278.75 

Au — Au — As 


6000.0 

Au — Ga — Ga 


1237.7 

Au — As — As 


5600.0 

Au — Ga — As 


3270.0 

Ga — Ga As 


1900.0 

Ga — As — As 


4600.0 
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Table 2. Calculated binding energies for varying numbers of gold 
atoms deposited on the (110) surface of GaAs. The total binding 
energy is denoted by E tot and e b represents the average binding energy 
per Au adatom. (Energies are given in eV). Letters a through g 
indicate for each cluster favorable positions of deposited Au atoms on 
the (110) surface as shown in Figure 1. 


Number of Sites 

Gold Atoms (See Figure 1) — Et 0 t 


1 

2 

3 

4 

5 


a 

a,b 

a,b,c 

a,c,d,e 

a,c,d,f,g 


2.69 

2.69 

5.44 

2.72 

7.44 

2.48 

10.40 

2.60 

12.90 

2.58 
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Figure 1 . A schematic top view of the (110) surface. Large and 
small open circles represent relative positions of Ga atoms located m 
the first and second layers. Large and small solid circles represent 
relative positions of As atoms in the first and second layers. Large 
open circles with letters are favorable sites for Au atoms of deposited 
clusters. A unit surface cell is indicated by dotted lines. 
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APPENDIX 2. 


A COMPARATIVE STUDY OF SILICON EMPIRICAL 
INTERATOMIC POTENTIALS 


1 


1 INTRODUCTION 


In recent years, there has been a surge of computer simulations for complex 
materials science-related phenomena, e.g., molecular dynamics simulations of 
melting, epitaxy, and crystal growth. In these types of calculations or 
simulations, it is, of course, desirable to use accurate first-principles quantum- 
mechanical methods. 1 However, because they require a large computational 
effort to accurately solve the Schrddinger equation, these methods are currently 
limited to studies of static properties for systems involving only a few tens of 
atoms. Nevertheless, progress, although slow, is being made to circumvent 
these limitations . 2 On the other hand, although they generally lack the accuracy 
of the former methods, empirical interatomic potentials can handle much larger 
systems and can be used to study static as well as dynamic properties of such 

systems. 

The theory of interatomic potentials for ionic systems and metals is rather 
well established as indicated by the remarkable success of some methods, e.g., 
the shell model 3 and the embedded atom method , 4 to accurately predict a wide 
range of properties. Unfortunately, the theory for covalent solids is less 
developed despite the many attempts of the last several years to model such 
strongly bonded materials. 

Not surprisingly, because of its technological importance, silicon has been the 
prototype material for developing empirical potentials. About sixteen such 
potentials have appeared in the literature in the last seven years. 5 * Of course, 
every new potential is claimed by its originators to be superior, i.e., more 
accurate and/or more transferable than its predecessors. While these claims 
are often valid to some extent, such improvements are almost always achieved 
by sacrificing other properties. 1CM3 Also, very often it is not truly clear what 
causes the better description. Is it due simply to a more flexible functional form 
and/or a better fitting strategy 10 or does the new potential really give a better 
description of covalent bonding? 12 Even the question of the range of interactions 
in covalent solids is not well understood and very often longer range forces 
are arbitrarily neglected for convenience. Due to the complexity of the structural 
chemistry of silicon and of the empirical nature of these potentials, the answer to 
these questions is certainly not an easy task. We attempt here to partly address, 
albeit indirectly, some of these issues by performing a comparative study of six 
of the aforementioned potentials. We will not consider here the so-called 
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valence-force potentials which can only describe small distortions from 
equilibrium. Stoneham, Torres, Masri, and Schober 25 have performed a 

comparison of eight such potentials for silicon. 

The potentials considered in this study are those of Pearson, Takai, 
Halicioglu, and Tiller (PTHT), 5 Stillinger and Weber (SW), 6 Biswas and 
Hamann (BH), 10 Tersoff (T2 and T3), 12 - 13 and Dodson (DOD). These 
potentials differ in their degrees of sophistication, functional form, and range of 
interactions; they thus constitute a good representative sample of existing 
potentials. Only a few comparative studies of some of these potentials have 
been performed. In general, these studies involved only a few of the potentials 
mentioned here and were limited to some specific properties of silicon. Khor and 
Das Sarma 26 used SW, DOD, and the first Tersoff potential (Tl) 8 in a study of 
Si(100) surface reconstruction. Li, Chen, Allen, and Broughton 27 investigated 
the energy and vibrational spectrum of the Si(lll) 7x7 surface using SW and 
T2. Halicioglu, Pamuk, and Erkoc 28 studied Si2 - Si4 and reviewed the results of 
some bulk and surface properties obtained with PTHT, SW, Tl, and DOD. 
B artel t, Williams, Phaneuf, Yang, and Das Sarmas 29 used SW, Tl, and T2 to 
model the orientational stability of silicon surfaces including (100), (111), (112), 
and (113). Finally, Bolding and Andersen 20 performed an extensive study of Si2 
- Si i o and compared the results obtained using their potential with those 
obtained via BH, SW, and T2. Their comparison for bulk and surface properties 

was more limited. 

Here, these potentials are thoroughly tested in a substantial region of 
configuration space, including small clusters, bulk, and flat surfaces. A great 
many results are new and, except for the bulk phonon frequencies, all those 
results already published by other researchers have been reproduced here for 
consistency and also because very often, they are incomplete. For instance, the 
calculations for the bulk monovacancy are more complete here since, in 
contrast to previous studies, we consider and report all possible configurations. 

The calculations for clusters (Si2 - Si6) were performed here with and 
without the cutoff function. Results obtained via PTHT and DOD for S 15 - Si6 
and via T3 for Si 3 - Si 6 are new. For Si 4 , the present work (with PTHT and 
DOD) is more complete than the study of Halicioglu, Pamuk, and Erkoc. 
Bolding and Andersen 20 performed calculations on Si2 - Siio using T3 but did 
not report any specific result. Moreover, although studies were previously 
performed using BH, SW, and T2, we identify here some new minima for S 14 - 
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Si6- We considered several bulk point defects: monovacancy and four types of 
interstitials, tetrahedral, hexagonal, bond-centered, and split or dumbbell. New 
results are presented for all these defect structures (PTHT and DOD), for the 
bond-centered interstitial (BH), and for the split interstitial (BH and T2). In 
addition to the three elastic constants of the cubic diamond structure, we also 
calculated here the pressure derivative of the bulk modulus, B’, and Kleinman s 
internal strain parameter, £. All calculations related to the elastic properties are 
new for PTHT, BH, and DOD. The results for B’ (SW, T2, and T3) and for C 
(SW and T3) are also new. For the bulk polytypes of silicon, we considered, in 
addition to the cubic diamond phase, all structures for which ab initio data is 
available. These are hexagonal diamond, BC-8, p-tin, simple hexagonal, the 
cubic phases, HCP, and the graphitic phase. In contrast to some previous 
studies, 30 the c/a ratios of the p-tin and all hexagonal phases were optimized. 
We also calculated the bulk modulus of those structures for which ab initio 
results are available, e.g., BC-8, p-tin, simple hexagonal, and the graphitic phase. 
The calculations performed with PTHT (hexagonal diamond, BC-8, and simple 
hexagonal), BH (graphitic phase), SW (BC-8, p-tin, simple hexagonal, and 
graphitic phase), DOD (BC-8 and simple hexagonal), T2 (BC-8), and T3 (p-tin 
and simple hexagonal) are new. The results for the pressure-induced phase 
transformations are new. Takai 30 reported that cubic diamond transforms under 
pressure to the P-tin phase but the c/a parameter of this structure was not 
optimized. Also, Biswas and Hamann 10 reported a transition under pressure 
from cubic diamond to simple cubic; we found here that this is only the second 
transition, the first involves a compressed HCP structure. All results for surfaces 
obtained with T3 are, in general, new. Other new calculations involving surfaces 
are: (110) (all potentials except SW); (100) 2x1 (PTHT); (100) c 2x2 (BH, 
DOD, and T2); (100) Pandey defect structure (PTHT, BH, and T2); (111) 2x1 
(SW); (111) 2x2 adatom covered structures (all potentials), (111) 
(2n+l)x(2n+l) DS and DAS structures with n = 1 - 4 (all potentials except for 
the 7x7 DAS surface with SW and T2). We also present results for the surface 
stresses of all these surfaces. In general, our results are in agreement with those 
already published by other researchers; however, there are some discrepancies, 
in particular with those involving BH. 

By a systematic comparison between these potentials, we identify similarities 
and differences and attempt to find and understand their ongins. We also assess 
the transferability and accuracy of these potentials with respect to experiment 
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and first-principles methods and discuss their limits and validity in quantitative 
modelling of materials phenomena. This paper does not, however, a ess e 
question of the theoretical justification or basis of interatomic potentials in 
semiconductors. This question has been dealt with by Carlsson. y 

highlighting the strengths and weaknesses of these potentials and by presenting 
such a large number of test results in a complete and clear manner, it is hope 
that it will help future users to select those potentials best suited for their needs 
as well as help future researchers desiring to either improve on these potentials 

or to develop new schemes. . , . M . 

The potentials are described in Sec. 2. The computational procedure 

discussed in Sec. 3. Results for small clusters, bulk phases, and flat surfaces are 
presented and discussed in Secs. 4, 5, and 6, respectively. Section 7 is a revtew 
of other potentials not considered in this work. Finally, a general discussion and 

conclusions are presented in Sec. 8. 

2 POTENTIALS 

Following Carlsson's classification.* PTHTBH.andSWwillt* retoedjo 

as cluster potentials and the last three as cluster functionals. SW, DOD, T2 and 
T3 are first-nearest neighbor models (second neighbor interactions are imp 1 y 
included in the bond bending term). BH and PTHT include interactions up to the 
third and seventh shell, respectively. We will sometimes refer to the fonner and 
latter group as the short-ranged and longer-ranged potentials, respective y. 

In order to have the units of energy and length in eV and A, respectively 
(and also for our own convenience), the notation used by the origin aut ors 
has been modified. It is, of course, straightforward to recover the ongina 
notation. Moreover, except for the PTHT and SW potentials the parameters 
listed in Table I are, or correspond to, those given in the original references, lb 
original parameters of PTHT and SW give a bulk cohesive energy for diamond 
silicon of -5.45 and -4.34 eV, respectively. It is convenient for comparison 
purposes to have a common basis which we choose as the expenmen al 
cohesive energy and lattice parameter of diamond silicon (the lat ‘ lce P arameter 
is predicted afmost exactly by all potentials). Thus for PTHT and SW bo h o 
the original two- and three-body energy parameters are multiplied by the smne 
scale factor of 0.85 and 1.07, respectively. The scaling of the energies in th s 
manner, does not affect the equilibrium structures previously determined. Note 
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that vibrational frequencies must be multiplied by the square root of these 
factors. 

2.1 Cluster Potentials 

The cluster potentials model bonding with classical two- and three-body 
potentials. The potential energy function gives the structural energy, E, which is 

written as, 


E = V 2 (r U 5 + X V 3 ( r ij' r ik - r ik > 

i,j Uk 

The primes indicate that all summation indices are distinct. Also because 
V 3 (i,j,k) is symmetric with respect to an interchange of i, j, and k, the more 
restrictive constraint, k > j > i, in the triple sum must be used with PTHT. The 
condition k > j is used with BH and SW because V 3 is symmetric in only j and k. 
The parameters listed in Table I correspond to these conditions. In general, the 

two-body potential is given by, 

V 2 (r) = f c (r)(A 1 (p 1 (r)- A 2 <p 2 ( r )) 

where f c is a cutoff function and <p s is a decaying function of r. The three-body 
potentials have different forms and will be given separately for each potential. 

2.1.1 The PTHT potential 

Pearson, Takai, Halicioglu, and Tiller 5 used the Lennard-Jones and Axilrod- 
Teller potentials for the two- and three-body terms, respectively. The potential is 
cut abruptly to zero at the cutoff radius, R c . In our notation, the functions are: 

f c (r) = 1 if r < R c 

= o otherwise. 

cp s (r) = r'^s s = 1,2 O) 


6 



V 3 (rij,rtk,rjk) = Z vfrij) y(r±) vj/(rjk) g(0i,0j,0k) 


V(r) 



g(0i,0j,0k) = 1+3 cos0i cos0j cos0 k 

0i is the angle subtended at atom i by atoms j and k; 0j and 0k are defined in a 

similar manner as shown in Fig. 1. 

There is, of course, no theoretical justification for using these potentials (in 
particular the Axilrod-Teller potential) to describe bonding in a covalent 
material. 5 These particular functional forms were chosen for purely pragmatic 
reasons and must, therefore, be viewed only as fitting functions. While the PTHT 
potential is not flexible enough (particularly with respect to the secon 
derivatives), it is appealing because it has only three adjustable parameters ( our 
if the cutoff radius is included; A-i and are fixed). This fact can be more 

appreciated when binary or ternary systems are considered. 

The parameters were fitted to a minimal database containing the on 
lengths of the dimer and trimer and the lattice parameter and cohesive energy of 
the'diamond structure. 30 The parameterization was done using infinite lattice 
sums (R c = ~). All results presented in this work correspond to R c - 7.3 A 
which yields results virtually indistinguishable from those obtained using an 
infinite cutoff radius. V 2 (r) is essentially zero at 7.3 A. Reducing R c to 5 A, i.e., 
including interactions up to third neighbors only, changes the lattice parameter 
by less than 0.5%, the bulk cohesive energy by 4% and the elastic constants by 
less than 10%. Similar changes are expected for surfaces, e.g., for the (111) lxl 
surface, the changes are 2%, 3%, and 9% for the first interlayer contraction, the 
surface energy, and surface stress, respectively. These changes are much 
larger if the potential is cut below the third shell. Also, the energetics of those 
surface defects which induce strong atomic distortions, e.g., surface vacancy 
and some types of steps, are much more sensitive to the cutoff radius. 

The PTHT potential has been extensively used to study bulk phase 
transitions, 33 surface reconstructions, 30 > 34 ’ 37 surface point defect formation and 
diffusion, 30 ’ 34 ’ 38 step reconstruction and interaction 3539 and their effects on two- 
dimensional nucleation. 39 The same potential has been extended to binary and 
ternary systems, e.g., GaAs, Si-GaAs, and Al-GaAs. 
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2.1.2 The Biswas-Hamann potential (BH) 


This potential is a simpler version of the original BH potential. The three- 
body potential is separable. Biswas and Hamann 10 used Gaussians for the radial 
functions and a Fermi-like function for the cutoff function. These functions are: 



f | 

f r-oY 

f c (r) = 

j 1 + exp 

1 1 1 


l 

l l 1 /J 

= 

0 



if r<Rc 
otherwise. 


<PsOO = e' Xs r2 


s= 1,2 


2 

V 3 (r ij ,r ik ,e i ) = X Z s m/ s (ry ) y s < r ik )8s ( e i > 


s=l 


¥s (r) = e - a s r2 f c (r) 


px 1 

g s (0i) = (cosGi - cos9 0 ) 

In the original paper, 10 a was written as R c suggesting that it is the cutoff 
radius; In fact, f c is 0.5 when r = o (for any value of \i) and 0.03 when r - 5.0 A. 
However, the results do not seem to be sensitive to R c as long as R c ^ o. For 
instance, the lattice parameter and cohesive energy of the diamond structure 
are 5.4126 A and -4.6267 eV when R c is equal to a; they are 5.4318 and -4.6045 
at 5.0 A and at 3.0 A they are 5.1565 and -4.9748. All the results presented 
here were obtained with R c = 5.0 A. V 2 (r) is essentially zero at this value. We 
have not used the separability of the three-body potential. Compared to the 
original notation, we have Z s = 2B S . The fitting database included the cohesive 
energies of a set of bulk phases, the formation energies of self-interstitials, and 
the surface energy of diamond Si(lll) lxl and of metallic simple hexagonal and 
cubic (100) surfaces obtained from an ab initio LAPW calculation. This 
potential has been used to study microclusters and bulk point defects, 
amorphous silicon, 43 cluster and atom deposition on Si(lll), 44 amorphous and 
epitaxial film growth on Si(lll), 45 and surface reconstruction. 3 
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2.1.3 The Stillinger- Weber potential (SW) 


The radial function, tp s . is given by (1). The other functions are: 


f c (r) = exp 


= 0 


'_k_' 

r-R c j 


if r < Rc 
otherwise. 


V 3 (rij,rik,rjk) = Z \|f(nj) g( 0 i> 


y(r) = [fc(r)] a 


g(0i) = (cos0i - cos0 o ) 2 

Compared to the original notation, 6 we have Ai = eABoP, A2 = eActf, Z = tK 
r c = ao and [l = g. Note that, besides acting as a cutoff function, f c defines the 
attractive branch of V 2 (since \ 2 = 0) and the radial functions of V 3 ; therefore, 
the results are very sensitive to variations in R c . Stillinger and Weber fitted the 
parameters to the lattice constant and cohesive energy of the diamond structure 
with the added constraint that the melting point and the structure of liquid silicon 
be well described. SW is by far the most widely used potential. It has been used 
to study clusters, 46147 lattice dynamics, 27t48 ’ 49 bulk point defects, the liquid 1 
and amorphous 52 ' 55 states, surface diffusion 56i5? and reconstructions, 
26.27,29,37.58.59 si(100) stepped surfaces, 60 the liquid- vapor 61 and crystal-melt 
interfaces, 51 ’ 62 pulsed melting of surfaces, 63 epitaxial growth from the vapor, 64 ’ 
66 liquid-phase epitaxy, 67 ' 69 and growth of amorphous films via atom deposition. 
70 This potential has been extended to Ge, 61>71 sulfur, 72 fluonne, and the Si-F 

system. 74 

2.2 Cluster Functionals 

All the cluster functionals considered here are of the Tersoff type. In this 
scheme, bonding is modeled with pairwise interactions but with the attractive 
term depending on the local environment which effectively includes many-body 
interactions. 8i ^ 2 The structural energy has the form. 
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( 2 ) 


E = f c (r 'j > ( A > <Pl (r ij ) ■ A 2 ^ (r *J } p( ^j ) ) 
« « 
l.J 


where p is a measure of the bond order and is a function of the effective 
coordination number, £ij, given by, 

CtrXvstm.ik.eo 

k*i,j 

V 3 ( r ij ’ r ik »®i ) = ¥( r ij » r ik )§(®i ) 

Note that, in general, p(Cij) * p(Cji)- The two-body energy can be extracted from 
E by rewriting (2) as 


E = -^'v 2 ( r ij) + ^X A 29 2(r ij )f c (r ij>( 1 - P(Cij)) 


(3) 


i.J 


i.J 


If [l - p(Cij)] is replaced with a Taylor expansion about some reference it can 
be seen how this scheme effectively includes many-body interactions. We will 

refer to the second term in (3) as the three-body energy. 

Brenner 76 showed that the Tersoff formalism is similar to the embedded 
atom method. 4 The two expressions for the structural energy can be made 
identical with the proper choice of functional forms and parameters. 


2.2.1 The Dodson (DOD) potential 

This potential is a simple modification of the first Tersoff potential. 8 The 
functions are: 


fc(r) 


_ 1 

( « 

fn(r-R c )Y| 


1 - cos 


” 2 

V 

l M- JJ 


= 1 


= 0 


if R c - |X < r < R c 

if r < R c - |i 
otherwise. 
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<ps(r) = e‘ X s- r 

p(Cij) = ex p['Cij n ] 


s = 1,2 


(5) 


¥ ( r ij * r ik ) = 


92 ( r ik ) f c( r ik) 
92( r ij) f c( r ij> 




(3 

^ “ -ScosGj 
T| + e 1 


To determine the parameters, Dodson used the lattice parameter and cohesive 
energy of the diamond cubic, simple cubic, BCC, FCC, and HCP structures, the 
bond length and energy of the dimer, and the bulk modulus of the diamond 
structure. To our knowledge, this potential is perhaps the least tested compared 
to the others. It has been used to study low-energy beam deposition of silicon 
and surface reconstruction. 9,26,37 

2.2.2 The Tersoff potentials (T2 and T3) 

T2 and T3 correspond to two different parameterizations of the same 
potential. f c and <p s are given by (4-5), respectively. The other functions are. 

p(Cij) = ( 1 + Cij n ) 2n 
¥(rij,rik) = fc( r ik)-exp[a 3 .(rij - rik) 3 ] 


g(0 i) P + g 2 g 2 + (cosOj — COS0Q 

The potential parameters were determined by fitting to a database containing 
the cohesive energy, lattice parameter, and bulk modulus of the diamond 
structure, and the cohesive energy of bulk polytypes of silicon. ’ For T3, an 
additional constraint was added in order to reproduce the three elastic constants 
to within 20%. 13 T2 and T3 were used to study microclusters, lattice 
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dynamics n - 15 - 27 bulk point defects, 12 - 13 - 78 the liquid 12 - 13 and amorphous 
states, surface reconstructions, 12.13,27,29,37 i ow . e nergy beam deposition on 
Si(100) 2x1. 79 This potential was also extended to carbon 80 and to the systems 
Si-C and Si-Ge. 8 '. 82 

2.3 Discussion 

The two-body functions V 2 (r) are plotted in Fig. 2. The open circles 
correspond to an accurate ab initio calculation (based on MRCI) of the energy 
of the ground-state of Si 2 . 83 The T2 and T3 curves are very similar for r < 2.7 
A. The attractive branches of the BH and SW curves are also similar over the 
range 2.4 < r < 3.0 A. While the shapes of these curves are different, in general, 
they all have about the same depth with the exception of DOD which is much 
stronger. The PTHT curve is the steepest with the largest curvature at the 
equilibrium bond length; this is reflected in the large values of the vibrational 
frequency of the dimer and the bulk elastic constants. Note the small bump on 
the curves of the cluster functionals; this is due to the abrupt cutoff function. 

The angular dependence, g(0), of the three-body potentials is shown m Fig. 3. 
The PTHT function is very different from the others; it is the only function which 
is negative for 0 > 117° resulting (in a few cases) in negative three-body 
energies for configurations with large bond angles. The DOD curve is a 
monotonic decreasing function of 0 with a minimum at 180° like the PTHT 
curve. The other curves have the same shape with a minimum at 0 = 6o- The 
T2 curve has a very shallow minimum at 0 = 90° and is symmetric with respect 
to that angle. For BH, since ai and a 2 are almost equal, gi(0)(l + (Z2/Zi)g 2 (0)) 
is plotted to show the effect of the cos 3 0 term which is negligible as shown in 
Fig. 3. A better comparison of the three-body potentials is provided in Fig. 4 
which shows the variations of the three-body energy with angle 0 for a triplet of 
atoms i, j, and k (Fig. 1) with ry = r ik = 2.351 A (the equilibrium bond length of 
the diamond phase) and 0i = 0. The three-body energy for the cluster functionals 
is defined as the second term of (3). What is actually plotted in Fig. 4 is the 
contribution of atom i to the total three-body energy of the triplet. There is no 
curve for T2 because this energy is essentially zero (on the scale of the plot) for 
the apex atom in an isoceles triangle; the total three-body energy does not, 
however, vanish (for 0 < 40°) because of the contribution of atoms j and k. One 
immediate consequence of this is that T2 will favor triplet configurations leading 
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to an equilateral triangle thus making the total energy almost completely 
controlled by the two-body potential. PTHT is the most repulsive for small 
angles less than about 45°, and SW and DOD are stronger than the other 
potentials for 0 > 135° and 50° < 6 < 140°, respectively. One very important 
result to be discussed later is that all these potentials handle rather well small 
angular distortions around the tetrahedral angle and, to a lesser extent, those 
leading to angles somewhat larger than 109°, but with the exception of T2, they 
completely fail when dealing with small angles (< 80°). This is fully reflected in 
the curves of Fig. 4. Finally, the large variations displayed by these functions 
makes a comparison between the potentials very difficult, indeed. 

3 COMPUTATIONAL PROCEDURE 

Unless indicated otherwise, all results presented here are static, i.e., T = 0 K. 
In the static limit, we use the so-called potential approximation where we only 
deal with the mechanical equivalent of the true thermodynamic properties. 84 
The total energy was minimized using a conjugate gradient technique. The 
minimization is stopped when the force on each atom is at most 0.001 eV/A 
(typically 10' 1 * * 4 eV/A). Note that, in general, energies converge more rapidly than 
stresses. For surface calculations, we used fixed atoms to simulate a rigid 
underlying substrate and used enough moving layers to minimize the interaction 
between the exposed surface and the moving-fixed interface. 

The surface energy per unit area, y, is defined as, 

y = 1(E-NE c ) (6) 

where E is the total energy for the simulation of N atoms with bulk cohesive 
energy Ec and A is the area of the exposed surface. 

The total stress T a p is given by 

1 3E 

IoP = va^ 

where V is the volume of the system, a and p represent any two of the 

cartesian coordinates and 8 0 p is the Lagrangian strain tensor. Since the total 
potential energy and the stresses are given as a sum of contributions of each 
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atom in the system, it is natural to define an atomic energy, e(i), and stress, 
XapO)> such that 

E = X^ and T ap = ^X* a P^ 

i i 

While the definition of the atomic energy and stresses appears to be arbitrary, 
the availability of the energy and stress distributions can be very helpful in 
analyzing complex defect structures. 34,82 Also, they are helpful in assessing any 
finite-size effect which can be present in the calculation. For example, well 
away from a defect, e.g., at a surface, the atoms must be representative of the 
bulk environment; this can be readily checked by looking at these distributions. 

A two-dimensional surface-stress tensor, G a p, can a ^ s0 ^ defined to further 
characterize the surface, 85,86 


CT ap ~ 


1 d(Ay) 

A 3e a p 


= y5 aP + 


3e ap 


(7) 


where 5 a p is the Kronecker delta. For a liquid, the strain derivative in (7) 
vanishes, and the surface stress, called surface tension, is numerically equal to 
the surface free energy. For a solid surface, d-y^ap can be positive or negative 
and can thus lead to a tensile or compressive surface stress. Also, as shown for 
the (100), (110), and (111) 2x1 surfaces, o a p is not necessarily isotropic. 

Unless indicated otherwise, all surface energies and stresses will be given in 
eV/lxl cell (for simplicity, we will sometimes drop the per lxl cell). To obtain y 
and G a p in eV/A 2 , divide by the area of the lxl primitive surface cell. Also, 
surface energy will refer to the value calculated from (6); the relative energy of 
a relaxed or reconstructed surface is, in general, its surface energy relative to 
that of the ideal lxl surface. 

For point defects, we calculated defect energies both at constant volume and 
constant zero pressure. For both cases, the reference is the perfect diamond 
lattice (pressure, P = 0; atomic volume, £2 = £2e) so that the formation energy, 

E f is given by, 


E/=E-(N± 1)E C 
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where E is the total energy for the simulation of (N ± 1) atoms, +1 for a single 
interstitial and -1 for a monovacancy. We used a large cell containing (980 ±1) 

atoms and the cartesian coordinate system was taken as x = [110], y-1001], 
and z = [110]. Periodic boundary conditions are applied in all three directions in 
order to simulate a bulk environment. For the constant volume calculation, the 
cell with a defect has to be compressed (vacancy) or expanded (interstitial) to 
bring the atomic volume to the equilibrium value, Qe- This is done by scaling the 
lattice parameter by the factor (1 ± 1/N) 1 / 3 . For N = 980, the change in lattice 
parameter is only ± 0.03% but it is ± 0.5% for N = 64. To simulate a constant 
pressure one would normally use molecular-dynamics or Monte Carlo methods. 
But these methods are much more time-consuming than a static method. If one 
is only interested in studying pressure effects at zero temperature on equilibrium 
structures, it is desirable to devise a static method where the volume is a 
variable. We have implemented the conventional method used in isobaric Monte 

Carlo simulation 87 into our static program code. 

In the "constant volume" static method, the objective function is the total 
energy E(r) and the variables are the 3N atomic coordinates (r) in the system 
of N atoms. The volume, V, of the computational cell is held fixed. In the 
"constant pressure" version, the enthalpy H{p,V) = E{p,V) + P ex V is minimized. 
There are now (3N + 1) variables which are V and the 3N dimensionless 
atomic coordinates, {p = V'^r). The externally applied pressure, P ex , is set to 
zero in this work. Note that, like the atomic coordinates, the periodic boundary 
vectors must be scaled by V- 1 / 3 at each iteration step. In practice, we found it 
necessary to constrain the volume to vary within a prescribed range. 

Because very often experimental data are lacking, people involved in 
atomistic computer simulation rely heavily on results obtained from first- 
principles calculations, not only for comparison purposes but also for use as 
input for the parameterization of the interatomic potential. In this work, only the 
properties of Si2 and those of the perfect diamond lattice are directly compared 
to experiment. Other bulk and surface properties are compared to results 
obtained from ab initio calculations using mostly the self consistent 
pseudopotential method within the local density functional theory (simply 
abbreviated DFT from now on). Results for Si3 - Si6 are compared with those 
obtained from an ab initio molecular orbital calculation. 88 
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Very often the ab-initio data are obtained from different sources, i.e., different 
techniques or calculations performed with different input calculational 
parameters. While the DFT method is certainly one of the most accurate 
theoretical tools for the calculation of the ground-state properties of solids, * it 
is important to keep in mind such differences in the data when making 
comparisons. Furthermore, for the parameterization of the interatomic potential, 
it is critical to use consistent and compatible data especially when used in 
conjunction with some experimental data, which is often the case. While several 
input calculational parameters can affect the result of the calculation (care is 
generally taken to minimize the effect of most of them) , the plane-wave cutoff 
energy, E pw , is the limiting factor. Structural parameters, e.g., the lattice 
parameters, are weakly dependent on E p w> but energies converge more slowly 
below about 10 Ry . 89 ' 92 Typically, E pw is about 6 Ry. The effect of E pw can be 
appreciated by considering the calculation of the energies of adatom-covered 
Si(lll) surfaces. Northrup 93 performed a calculation at 6 Ry and he 
successfully determined that the T 4 adsorption site is energetically favored over 
the H 3 site (both sites have three-fold symmetry but the T 4 adatom has a 
second-layer atom directly below in contrast to the H 3 adatom which does not). 
He also found that the ^3x^3 structure was more stable than the 2x2 structure. 
In fact, the opposite was found to be true by Meade and Vanderbilt who 
performed a more accurate calculation at 12 Ry . 91 Other examples are provided 
by the formation energies of bulk point defects in silicon for which extensive 
calculations lead to uncertainties of up to 2 eV 21 > 32 ’ 94 ' 96 and by the energy of 
buckled and symmetric dimers on the Si(100) surface for which seemingly 
similar calculations lead to opposite results (see Sec. 6 ). 

4 CLUSTERS 

The only accurate experimental data for silicon clusters are the bond length, 
bond energy, and vibrational frequency of Si 2 92 and the binding energy of Si3 • 
No experimental data on the structures of Sin (n ^ 3) is available. However, 
accurate ab initio calculations have been performed on microclusters. 88 ’ 98 
Raghavachari (KR) performed an ab initio molecular-orbital calculation and 
considered several configurations and electronic states for each Si n cluster (n = 
2 - 6 ) in order to determine the ground-state structures . 88 The ground-state 
structures were found to be somewhat compact with an average bond length of 
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about 2.28 A (shorter than the equilibrium bond length of the diamond structure) 
and bond angles in the range 60° - 90°. These structures are significantly 
different from the configurations derived from microcrystal fragments. These 
latter lie much higher in energy than the former. The lowest-energy structures 
for Si 3 - Si 6 are the isoceles triangle, rhombus, trigonal bipyramid, and edge- 
capped trigonal bipyramid, respectively. Each cluster can be formed from the 
previous one by attaching an extra atom at an edge- or face-capped site. These 
Sin structures are three-dimensional in nature for n > 5. At the level of theory 
considered and based on the experimental data of Si2 and Si 3 , KR estimated that 
only about 80 % of the binding energy of these clusters is recovered. He, thus, 
scaled the binding energy of all the clusters by the same factor of 1.2. The 
energies given in Tables II-III are the scaled values. 

Several studies on small silicon clusters were performed using some of the 
potentials considered here. Halicioglu, Pamok, and Erkoc 28 used the PTHT, SW, 
and DOD potentials for a limited study on Si2 - SU- Biswas and Hamann 10 used 
their potential along with a combination of steepest descents and simulated 
annealing techniques to determine the low-energy structures for Sin clusters (n 
= 3-6, 10, 32). A molecular dynamics simulation of Sin clusters (n = 3-17, 32) 
was performed by Blaisten-Baroja and Levesque using the SW potential . 46 Both 
neutral and positively charged clusters were examined. Feuston, Kalia, and 
Vashista 47 also used the SW potential to study the fragmentation of Si n clusters 
(n = 2-14) with a molecular dynamics technique. Bolding and Andersen 
performed a calculation on Si2 * Si i o using BH, SW, T2, and their own potential. 

We have performed our own calculations on Si2 - Si6 using all six potentials 
along with a static method. For each potential, except SW, calculations were 
carried out with and without the cutoff function. Because simple minimization 
energy techniques are not adequate for finding global minima, we have 
considered many configurations for each cluster; this includes all geometries 
considered in previous studies as well as several other structures such as a 
planar C2V form for Sis (edge-capped rhombus) and a Csv form for Si6 
(pentagonal pyramid) as illustrated in Fig. 5. Moreover, many asymmetric 
structures derived from each of these configurations were optimized. Our results 
of the calculations with the cutoff function are summarized in Tables II-III. 
Since the energy of atoms infinitely separated is taken as zero, the binding 
energy, Eb, is the absolute value of the total potential energy as given in Sec. 2. 

4.1 Si2 and Si 3 
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The equilibrium bond length, r e , binding energy, D e , and vibrational frequency, 
0) o of the Si2 dimer, as obtained from the six potentials are shown in Table II. 
Most potentials used r e in the fitting database; it is, therefore, well reproduced 
by most of them. The largest discrepancy of 5% in r e occurs with SW. D e is less 
well described with a discrepancy ranging from 12% for DOD to -29% for SW. 
The best description of co 0 is obtained with DOD; PTHT overestimates it by 
50% and the other potentials underestimate it by about 10%. 

The ab initio calculations predict that the lowest-energy configuration for S13 
is an isoceles triangle with an apex bond angle of 77.8° and bond length of 2.165 
A which indicates strong multiple bonding character. 88 It was also shown that 
the linear structure is only a saddle point on the potential energy surface of S13 
and that there is also a low-lying state which lies only a few kcal/mol higher in 
energy; its structure is an equilateral triangle with bond length of 2.263 A. None 
of the potentials predicts the correct ground-state structure of Si3 (see Table II). 
As expected from Fig. 4, PTHT, T3, and DOD predict a miner with bond angle 
of 180°, 126.75°, and 180°, respectively, i.e., the angle at which the three-body 
energy is a minimum. These three potentials correctly predict the second 
minimum; they, however, overestimate the bond length. Not surprisingly, as 
noted in the previous section, T2 predicts an equilateral triangle as the ground- 
state structure of Si3. SW and BH were expected to predict an isoceles triangle 
with a bond angle equal to the tetrahedral angle as the lowest-energy structure 
which turned out to be an equilateral triangle. BH like T2, does not, however, 
predict a second minimum. With T2, the bond length and binding energy of Si 3 
are constant at 2.31 A and 5.25 eV for angles larger than 88°. The calculation 
performed without the cutoff function shows that, for this potential, there is, in 
addition to the same global minimum, a shallow local minimum corresponding to 
a linear structure with Eb = 5.2463 eV and a bond length of 2.31 A (the 
maximum between these two minima occurs at 120° with Eb = 5.2419 eV). 

4.2 Si4, Sis, and Si6 

The first three lowest-energy structures for Si4 ■ Si6 are presented in Table 
III. Some of the structures are illustrated in Fig. 5. The hexagonal chair and the 
chain are crystal fragments when 6 = 109.5°. In the corner-capped triangle, r34 
is slightly different from ri3. The rhombus in the corner- and edge-capped 
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rhombus is somewhat distorted. For both the trigonal bipyramid and square 
pyramid, there are two different geometrical arrangements, a flattened and an 
elongated form. ri 5 < m in the flat square pyramid and ris > H2 in its 
elongated form; atoms 1, 2, and 3 are bonded to each other in the elongated 
trigonal bipyramid while they are not in the corresponding flattened form. The 
wedge is a symmetrical stacking of two equilateral triangles. The orthorhombic 
bipyramid can be viewed as two edge-sharing distorted tetrahedrons. With 
PTHT, atoms 1 and 2 are bonded to each other in the asymmetric structure. For 
the linear structures, the numbering of atoms is from one end to the other as in 
Ref. 88. 

Let us first discuss the ground-state structures. Only PTHT correctly 
describes the ground-state structure of SU, i e., the rhombus. This structure is 
predicted to be a local minimum on the potential energy surface by DOD and 
T2. With BH, SW, and T3, it is not a minimum. For Sis, the quantum ground- 
state structure, e.g., the flat trigonal bipyramid, is a local minimum with all 
potentials (E B = 10.78 eV with PTHT and 12.03 eV with DOD). The edge- 
capped trigonal bipyramid, which is the global minimum on the quantum 
potential energy surface is only a local minimum on the surface generated with 
PTHT (Eb = 14.18 eV), BH (15.26), SW (14.95), DOD (14.89), T2 (23.08), and 
T3 (13.13). For the first three potentials, the structural parameters of the 
optimized structures are very different from those of the quantum structure. 

For a potential to be useful in studies of clusters, it should, as pointed out by 
Bolding and Andersen, 20 give at least a fair representation of the entire potential 
energy surface. That is, it should not only fairly descibe the energies and 
structures of global and local minima but also, and perhaps more importantly, not 
predict spurious minima (minima which do not exist on the quantum potential 
energy surface). In addition to the global minima, there are three, three, and one 
known local minima on the ab initio potential energy surface for SU, Sis, and SU, 
respectively. Not listed in Table III are the pyramid (Eb = 9.14 eV) for SU, and 
the pentagon (Eb = 12.10 eV, r = 2.39 A) for Sis. KR also reported that, for Sis, 
the linear structure and the tetrahedral crystal fragment lie higher in energy than 
all the other structures and that the former has a larger Eb than the latter. 88 It is 
not known whether these two structures are minima. 

For SU, all potentials predict that the tetrahedron is a minimum while the 
D2d structure is not. The pyramid is predicted to be a local minimum only by 
SW; however, its apex angle is 109.5° (a crystal fragment) compared to the 
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value of 78° for the quantum structure. This structure is degenerate in energy 
with the chain (another crystal fragment). For Sis, structures that are predicted 
to be minima are: the elongated trigonal bipyramid by PTHT (Eb = 10.75 eV), 
T2, and T3 (9.83 eV); the flat square pyramid by T2; and the pentagon by all 
potentials (for T2, E B = 13.12 eV, r = 2.31 A). For Si 6 , only BH and SW predict 
that the hexagonal chair is a local minimum; however, the bond angle is 107.2 
and 109.5°, respectively. The value for the optimized quantum structure is 93.6°. 

We now consider the spurious minima. KR determined that the following are 
not minima on the quantum potential energy surface: for S4, the comer-capped 
triangle (E B = 10.63 eV), the square (10.61 eV; r = 2.32 A), and the linear 
structure (8.75 eV); for Sis, the elongated square pyramid (16.16 eV; m = 2.30 
A, ris = 2.50 A); for Si6, the face-capped trigonal bipyramid (21.87 eV) and the 
tetragonal bipyramid (21.48 eV). The following structures are thus spurious 
minima: the square, the corner-capped triangle, and the linear structure (all 
potentials); the elongated square pyramid (all potentials except T2); face- 
capped trigonal bipyramid (DOD, T2, and T3); the tetragonal bipyramid (T3). 
Finally, we should mention that, with all six potentials, we have found many 
more local minima (about 15 overall for each potential) which, in general, are 
close in energy. 

PTHT and DOD favor planar structures. In fact, PTHT and DOD give rather 
similar descriptions of the structures of the small Sin clusters. This is consistent 
with the similar monotonic angular variations of the three-body energy for 
relatively larger angles. The most similar potentials are BH and SW. Both 
predict the same ground-state structure for each cluster. BH gives slightly larger 
binding energies and smaller bond lengths than SW. T3 is also close to BH and 
SW. It predicts the same ground-state structures (with similar bond lengths and 
binding energies) for SU and Sis* There are more similarities between T3, BH, 
and SW if all minima for Si4, Sis, and Si6 are considered. That BH and SW and, 
to lesser extent, T3 are similar can also be traced back to the behavior of the 
corresponding three-body potentials as shown in Fig. 4. Note the tendency for 
T3 to also favor planar structures. As expected from our earlier discussion, T2 
favors structures which lead to triplets forming an equilateral triangle or close to 
it. For example, for S4, the tetrahedron and rhombus are actually 4 and 2 edge- 
sharing equilateral triangles, respectively. All three minimum energy structures 
and also the trimer structure have the same bond length of 2.31 A, i.e., the 
equilibrium bond length of the dimer. It is this discrimination in favor of 
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structures having equilateral triangles as their building blocks which leads to 
three-dimensional configurations in agreement with the ab initio calculations. In 
fact, T2 is the potential that provides the fairest overall agreement with the ab 
initio calculations. 

In order to compare the cluster binding energies, a plot of the binding energy 
per atom versus the number of atoms in the cluster is shown in Fig. 6. The 
curve labeled KR corresponds to the scaled energies of the ab initio calculations. 
The similarity noted above for BH, SW, and T3 is also apparent in the binding 
energy. The PTHT curve is very similar to that of these potentials. The DOD 
potential, which showed similarities with PTHT for the structures of the Si n 
clusters, leads to larger binding energies than PTHT because its two-body 
potential is stronger. Compared to the scaled ab initio results, PTHT, BH, SW, 
DOD and T3 generally underestimate the binding energy of the Sin clusters 
while T2 overestimates it. Note that, for T2, the energy per atom is increasing 
with a relatively large slope and it is already 4.42 eV for n=6, i.e., very close to 
the bulk cohesive energy. This is a direct consequence of the fact that the bond 
bending forces are very small; thus T2 favors close packed structures because 
the total energy is almost totally controlled by the two-body potential. A best fit 
with the KR curve can be obtained for the binding energies by using a scaling 
factor of 1.43, 1.40, 1.47, 1.15, 0.85, and 1.39 for PTHT, BH, SW, DOD, T2, and 
T3, respectively. However, this will also change the bulk and surface energies 
and it will not change the fact that the equilibrium configurations are, in general, 
in disagreement with the ab initio results. To compare the relative stability of 
these clusters, and thus look at the possibility of magic numbers, one needs to 
perform a calculation of the fragmentation energy, E f r . This is the smallest 
energy involved in the dissociation of Sin into Sin-m + Sim- An accurate 
determination of Efr involves the investigation of all possible fragmentation 
channels. KR determined that E fr corresponds to the process Si n -> Si n -1 + Si 
and confirmed the presence of the magic numbers 4 and 6. 88 Using SW, 
Feuston, Kalia, and Vashista 47 also found the same fragmentation process for 
Si n (n = 2 - 14) and that SW do give the magic numbers 4 , 6, and 10. The kind 
of extensive study performed in Ref. 47 is beyond the scope of this work. 
However, assuming the same process as in Refs. 47 and 88, we find, for Si2 - 
Si 6 , that Si 4 is more stable than the other clusters with BH, SW, and T2. PTHT 
and T3 give more stability to Si 5 . With DOD, no cluster shows extra stability as 
indicated by the almost linear curve in Fig. 6. 
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All the results presented thus far were obtained with the cutoff function 
included with the potentials. We have also performed the same calculations with 
no cutoff; the cutoff function in SW is an integral part of the potential. The two- 
body potential function without the cutoff function , i.e., q(r) = V 2 (r)/f c (r), is 
essentially zero at r = R c for PTHT and BH. q(R c )/q(re) is 0.42, 0.47, and 0.54 
for DOD, T2, and T3, respectively. As expected, there is very little or no change 
in the results obtained with PTHT and BH. For the cluster functionals, there are 
little changes in the structural parameters, somewhat larger variations in the 
binding energies (particularly for T3), and some minor changes in the relative 
position of the minima. The global minimum of SU is now the rhombus for T3 
and that of Sis is now the flat square pyramid for T2. The largest variations 
were obtained with T3 but they are not large enough to change the overall 
picture presented thus far. 

Finally, It should be mentioned that, while our results agree with the works o 
Refs. 20, 28, 46, and 47, there is some conflict with the work of Biswas and 
Hamann. 10 First, the binding energies listed in Ref. 10 are consistently slightly 
larger (less than about 0.1 eV/atom) than the values presented here. The 
structure of Si 3 is in total disagreement; they found an isoceles triangle with 0 - 
79° r = 2.29 A, and Eb = 5.10 eV. Our result agrees with that of Ref. 20. With 
BH,’ we found all minima reported in Refs. 10 and 20 and many more. Our 
results for T2 are in agreement with those of Ref. 20 and here again we found 

other minima. 

In summary, the ab initio calculation 88 predicts that, for Si3 - Si6, there are 
overall twelve structures which are minima and seven which are not (there are 
probably many more than that). Considering these nineteen structures and using 
a simpler version of the rating scheme of Bolding and Andersen, 20 a potential 
will predict (m,n) minima with 0 < m < 12 and 0 < n < 7 being the number of 
correct and spurious minima, respectively. The rating of the ab mitio calculation 
is (12,0). It is (8,4), (8,4), (6,4), (7,5), (7,6), and (6,6) for T2, SW, BH, PTHT, T3, 
and DOD, respectively. The most serious limitation of these potentials is that 
they predict many spurious minima which are either global or close in energy to 
the correct global and local minima. A positive note is that, in general, these 
potentials do predict, like the ab initio calculations, that the structures derived 
from crystal fragments are not energetically favorable even though the 
potentials were built from crystal data. Within the framework of classical 
interatomic potentials applied to covalently bonded materials, it is the delicate 
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balance between the radial and bond bending forces which determine the 
equilibrium structure and the energetics of any system of atoms, e.g., clusters, 
bulk phases, or surfaces. The fact that these potentials do not describe correctly 
the equilibrium structures of the silicon clusters indicates that such a balance is 
not adequate at this point. 


5 BULK PHASES 


5.1 Crystal Stability 

In addition to the cubic diamond structure (a = 5.429 A; Ec = -4.63 eV/atom), 
silicon may exist in several simple and complex metastable structures. 89 These 
phases have been observed experimentally and most of them result from 
pressure-induced phase transformations. They are: hexagonal diamond which 
has the same density as cubic diamond (c/a = 1.653 ; a = 3.80 A), p-tin (c/a - 
0 552 * a = 4 686 A), BC-8 (x = 0.1003 ; a = 6.636 A), simple hexagonal (c/a = 
0.94 ; a = 2.527 A), HCP (c/a = 1 .698 ; a = 2.444 A), 99 and FCC. 100 Hexagonal 
diamond is formed with a combination of high-pressure and heat treatments. 
The BC-8 structure is observed upon unloading to atmospheric pressure from 
the high-pressure P-tin phase. 101 The other phases are the result of high- 
pressure phase transformations in the pressure range 0-800 kbar. - They 
occur in the sequence cubic diamond ~> p-tin -> simple hexagonal ~> HCP - 
> FCC with increasing pressure. 

Accurate DFT calculations have been performed on these structures as well 
as some other hypothetical phases, e.g., simple cubic, BCC, and graphitic 
structure. 99 The DFT database has been very useful not only as input for the 
parametrization of the potentials but also for comparison purposes when 
experimental data are unavailable as is often the case. To test for crystal 
stability, and as further comparison between the potentials, we performed 
calculations for all the structures mentioned above as well as for several two- 
dimensional structures. The axial ratio, c/a, of hexagonal diamond, p-tin, simple 
hexagonal, HCP, and graphitic silicon as well as the internal parameter, x, of 
BC-8 were optimized. The results for the optimized structures are shown in 
Table IV. In the DFT calculations the c/a ratios of the hexagonal diamond, HCP, 
and graphitic structures were not optimized. 89,102 
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All potentials, but PTHT, predict the cubic diamond phase as the most stable 
structure. The major result here is the unfortunate finding that with PTHT the 
lowest-energy phase is the simple hexagonal structure instead of cubic diamond. 
The axial ratio of 2.87 is so large that there is negligible interaction between the 
hexagonal layers. The second lowest energy structure is not even cubic 
diamond but a squared two-dimensional structure with a = 2.32 A and AE — Ec - 
E c (cubic diamond) = -0.22 eV/atom. This finding stresses the need, when testing 
for crystal stability, for considering all plausible phases including planar 
structures. When PTHT was first developed, crystal stability was tested with a 
minimum number of structures. 30 This situation is certainly not unique. As an 
example, the same problem occurred with the first potential developed by 
Tersoff. 8 Despite this pathology, we will continue in the next sections to present 
results obtained with PTHT because, as we will show throughout, this potential 
yields in most cases similar results to those obtained with DOD. 

Crystal stability is the first requirement any potential must fulfill for it to be 
useful. 104 If the potential gives as the most stable structure a phase other than 
the experimentally observed one, it cannot be used, in general, and particularly 
for melting and growth simulations. However, it can, perhaps, be used for a 
limited number of structural calculations in regions of phase space away from 
the pathological configuration; but one has to remain skeptical of such 
calculations. When the most stable structures are planar structures, as in the 
case with PTHT, one can easily foresee situations where the use of such a 
potential could lead to trouble. For example, in an unconstrained simulation 
involving (100) or (111) surfaces and where the lateral dimensions are relaxed, 
the layers making up the slab could break away if the temperature is high 
enough to overcome the energy barrier preventing bond breaking. In fact, this 
pathology is perhaps responsible for the results of the simulation of thin 
amorphous silicon films on crystalline silicon substrates performed by Erkoc, 
Halicioglu, and Tiller 105 They found that the dominant structural feature was a 
dense free surface skin with a void layer underneath for Si (100) and (111) 
substrates. 

A good fit to the energy of the hexagonal diamond phase is important 
because there are direct implications on the energies of stacking faults and of 
the Si (111) 7x7 surface. The DFT calculation shows that AE = 0.016 eV/atom 
in agreement with the observation that the cubic and hexagonal diamond phases 
are closely related and with the small experimental stacking fault energies in 
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silicon along the <11 1> direction, e.g., 50 - 60 erg/cm 2 . 106 BH gives the best fit 
with a c/a ratio slightly greater than the ideal value in agreement with the 
experimental value. PTHT gives a vanishingly small AE. The other short-ranged 
potentials, being first-nearest-neighbor models, give, as expected, a zero AE. 

The BC-8 structure is of particular interest because it provides information 
about bond-bending forces. This phase has a BCC structure with 8 atoms per 
unit cell; it has two structural parameters, the cubic lattice parameter, a, and the 
internal parameter, x. It consists of distorted tetrahedra with small changes in 
bond lengths. As in cubic diamond, each atom in BC-8 has four neighbors but 
there are two different types of bonds resulting in two slightly different^bond 
lengths and thus two distorted bond angles, about 99° and 118°. ’ All 

potentials describe the structural parameters fairly well; however, with the 
exception of T2, they overestimate the energy. T2 gives a very small relative 
energy implying very weak bond-bending forces as already indicated in the 
previous section. SW, DOD, T2, and T3 correctly predict the bulk modulus; 
PTHT and BH overestimate it just as they do for cubic diamond (see Table V). 

The (3-tin structure is also of interest because it is the first phase the diamond 
structure transforms to under pressure. This phase has four nearest-neighbors 
and two second-nearest neighbors at a slightly larger distance making the 
effective coordination number, Z e ff, equal to 6. SW gives the best overall 
description (excluding the bulk modulus). Only T3 correctly predicts the bulk 
modulus. All but SW, underestimate the c/a ratio and only BH and SW describe 


the energy fairly well. 

The graphitic phase is important because (i) it is the only undercoordinated 
structure (with respect to cubic diamond) and (ii) it gives a measure of both the 
tendency of silicon for rehybridization from sp3 to sp 2 as observed on some 
surfaces and also of rc-bonding which was found by first-principles calculations 
to be weak in silicon compared to carbon. 102 Note again that in the DFT 
calculations, the c/a ratio was not optimized. 89 The value of 2.726 for graph, te 
was used. Our results for this value are AE = 0.34, 1.09, 1.27, 0.39, 0.70, and 
0.72 eV/atom and a = 3.92, 3.99, 4.10, 3.95, 4.01, and 3.99 A for PTHT,BH, SW, 
DOD T2 and T3, respectively. Note the excellent agreement obtained with T2 
and T3. The optimized structures, as obtained with these potentials have a much 
smaller axial ratio of about 1.2, remarkably about the same for all potentials. The 
energy is also comparable for PTHT, BH, and DOD and for T2 and T3. Zeff is 5 
for the atom having neighbors directly above and below in the adjacent layers. 
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PTHT and BH predict a compressed HCP structure with an axial ratio of 
0.59 and 0.69, respectively. In this structure, Zeff is 8 instead of 12 for the ideal 
HCP structure; this structure is very close in energy to the cubic and hexagona 
diamond phases. SW also predicts a HCP structure with a c/a ratio of 0.884 
smaller than the ideal value of 1.633; Zeff is only 6 in this case. For the ideal 
HCP structure the results are a = 2.81, 2.88, and 2.93 A and AE 0. , • . 

and 0.42 eV/atom for PTHT, BH, and SW, respectively. The lower energy of the 
compressed structure is due mainly, particularly for PTHT, to a lower three- 

The energy of all the structures is compared to the DFT results in Fig. 7. In 
this figure, the energy of the graphitic structure corresponds to the non 
optimized c/a ratio of 2.726. Also, for PTHT, BH, and SW, the energy of HCP 
corresponds to the ideal c/a ratio. Only DOD shows the same trend m energy 
(up to the HCP phase) as the DFT results. This is not surprising because 
Dodson included more structures in the fitting database than the other potentials. 
PTHT does a poor job in describing the energy of these phases. Only SW, D 
T2 and T3 predict the BC-8 structure as the next-lying phase after cubic and 
hexagonal diamond, and only SW and DOD predict p-tin as the fourth phase. All 
potentials correctly predict that the equilibrium atomic volume of cubic diamon 
is larger than that of the other phases (excluding the non-optimized graphitic 
structure). They also predict the increase of bond length with increasing 
coordination (the dependence is approximately logarithmic). However, the bon 
lengths are, in general, somewhat larger (particularly with SW and BH) 
compared to the DFT results. 

5.2 Phase Transformations 

As a further test and comparison, we finally discuss pressure-induced phase 
transformations. We have studied all possible phase transformations (mainly 
from the cubic diamond structure to all the other bulk phases) using t e 
procedure outlined in Ref. 89. Only T3 predicts correctly that cubic diamond will 
first transform to the p-tin phase at a transition pressure of 127 kbar. The 
transition volumes (normalized to the experimental equilibrium volume of cubic 
diamond) are 0.903 (cubic diamond) and 0.715 (p-tin). This agrees fairly we 
with the experimental values of 88-125 kbar, 108 0.918, and 0.710, respectively. 
T3 also predicts a p-tin to BC-8 transformation at 47 kbar compared to the D 
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result of 76 kbar. 101 The first transformation from cubic diamond is predicted to 
be to the compressed HCP phase by PTHT and BH, to BC-8 by SW and T2, 
and to simple hexagonal by DOD. The cubic diamond to fi-tin phase 
transformation is predicted to be (in the order of increasing transition pressure) 
the third transformation by PTHT (286 kbar, 0.93 (cubic diamond), 0.77 ((3-tin)), 
the fifth by BH (155, 0.93, 0.84), the second by SW (217, 0.86, 0.82), the fourth 
by DOD (205, 0.85, 0.76), and the fifth by T2 (270, 0.84, 0.73). Biswas and 
Hamann 10 found that cubic diamond would make a transition under pressure 
first to simple cubic. We found this transition to occur after the cubic diamond to 
HCP transition using their potential. 

5.3 Elastic and Vibrational Properties 

The description of the elastic properties (in particular, the shear constant 
C 44 ) constitutes a stringent test for the potentials. We have calculated the 
elastic constants using the homogeneous deformation method. 84 The results are 
presented in Table V. The cluster functionals were all fitted to the bulk modulus, 
B; it is thus well described. They also descibe well the pressure derivative of the 
bulk modulus, B'. SW provides a good fit to B but significantly underestimates B . 
PTHT and BH overestimate C 11 , C 12 , B', and B. DOD and T2 underestimate 
Cn by about 28%; SW and T2 overestimate C 12 by about the same amount. All 
potentials underestimate C 44 ; the worst fit is provided by T2, PTHT, and BH. 
For T2, this indicates once again the very weak bond-bending forces. Note that 
they also underestimate the second shear constant (C 11 - C 12 ). Also, none of 
the potentials correctly describe the negative Cauchy discrepancy (C 12 - C 44 ). 
PTHT and DOD overestimate significantly Kleinman’s internal strain parameter, 
£. 23 This is reflected in the large value of (C°44 - C 44 X C°44 is the value of C 44 
in the absence of internal displacement. 110 BH, SW, and T3 provide a good fit to 
C °44 and This shows that such a good fit does not necessarily result in an 
accurate value of C 44 . Tersoff developed T3 in order to improve on the elastic 
constants. T3 indeed describes the elastic constants better than does T2. SW 
also gives a good description. This is in fact quite remarkable considering that 
SW was not directly fitted to any elastic constants. For all potentials, the elastic 
constants satisfy the mechanical stability conditions indicating that the^cubic 
diamond structure is stable against all elastic homogeneous deformations. 
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The vibrational properties influence small distortions from equilibrium and are 
thus also important. In Table V are listed the phonon frequencies corresponding 
to four modes, the transverse acoustic at X, TA(X), the transverse optic at X, 
TO(X), the longitudinal optic and acoustic at X, LOA(X), and the longitudinal- 
transverse optic at T. The phonon spectrum for DOD is not available. BH gives 
the best overall description of these phonon frequencies. 10 For PTHT, despite a 
very poor description of the elastic constants, the phonon frequencies agree 
fairly well with experiment with the exception of vto(X). Note the almost 
perfect agreement for v TA (X). The full phonon spectrum has been determined 
by Pearson. 34 At the zone boundary, the acoustic modes are well described. 
However, the optical modes are too stiff; besides a larger zone center 
frequency, the modes increase with increasing wave vector, instead of 
decreasing. The fact that PTHT and BH give a better agreement with 
experiment for the transverse acoustic mode at X than the other shorter-range 
potentials is consistent with the observation that long-range interactions are 
necessary for an adequate description of this mode. 113 While SW also does a 
fairly good job with these phonon frequencies, the spectrum extends to higher 
frequency compared to experiment. T2 gives a poorer overall description of the 
phonon spectrum than SW. v T a(X) is underestimated significantly. T3 gives a 
description comparable to SW. 13 

5.4 The Universal Energy Relation 

All properties of the cubic diamond structure presented thus far probe only a 
very small region around equilibrium. It is desirable to seek information about 
the behavior of the potentials in regions well away from equilibrium. Recently, 
Rose, Smith, Guinea, and Ferrante 114 showed that the binding energy versus 
distance relation for a condensed system, independent of its bonding character, 
could be described by a universal energy relation (UER) given by, 

E(a*) = Ec (1 + a* - f 3 a* 3 ) e‘ a * 

Where a* = -n(r/r e - 1) and Tl = (9BQ e /IEc0 1/2 . Here r is the first nearest- 
neighbor distance and Q is the atomic volume while the subscript e indicates 
equilibrium values. The coefficient f 3 was fixed to a value of 0.05 which was 
determined from the thermal expansion of copper. 114 This led to a value of 4.74 
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for B’ in silicon compared to 4.2 for the experimental value. A value of 0.0052 
for f 3 leads to an exact fit whereas f 3 = 0 gives B' = 4.25; f 3 was set to zero in 
this work. In general, good agreement was found between the universal energy 
relation and experimental data as well as the first-principles calculations. The 
UER not only accurately models the region of moderate uniform compression 
but it also provides information about the behavior of the crystal under moderate 
uniform expansion where experimental data is lacking. 

The energy versus distance curves are shown in Fig. 8. Compared to the 
UER curve, BH and PTHT indeed do a poor job. The repulsive branch of the 
curve is well described by the short-range potentials. The cluster functionals 
also do a good job for expansion up to a bond length of about 2.8 A. Note that 
the UER suggests that the potential is long-range; however, the accuracy of the 
UER has not been demonstrated for large expansions. According to first- 
principles calculations, 83 the range of the two-body potential is about 5A (see 
Fig. 2). Moreover, the many-body term should fall off much more rapidly with 
distance than the two-body term. It thus seems reasonable to assume that the 
range of the potential should be no more than about 5 A. 

To further characterize the behavior of the energy-distance curve, we 
consider a theoretical property defined in Ref. 114 as a limit on the tension at 
which the crystal would rupture. This negative pressure, Pr, is the value at the 
minimum of the pressure versus distance equation. It was found that, for most 
crystals, Pr is typically 10 to 20% of the bulk modulus. Compared to the UER 
value of -0.16 Mbar, P R is -0.17, -0.14, -0.16, -0.34, -0.36, and -0.53 Mbar for 
PTHT, BH, SW, DOD, T2, and T3, respectively. The corresponding strain, e R , is 
9, 13, 20, 26, 26, and 20%, respectively while the UER value is 15%. The large 
overestimate of Pr and e R exhibited by the cluster functionals is due entirely to 
the short range and abrupt cutoff (Figs. 2 and 8). 

5.5 Bulk Point Defects 

In general, point defects involve large atomic displacements and rebonding 
around them. They provide information about bond-breaking energies and are 
important for testing the ability of the potentials to model such large atomic 
displacements. Several calculations or simulations of intrinsic defects were 
previously performed using BH, SW, T2, and T3. Using SW, Batra, Abraham, 
and Ciraci 50 performed an extensive molecular-dynamics simulation (at zero 
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pressure and temperature and in a cell containing 800 atoms) for four types of 
self-interstitials: the tetrahedral (It), the hexagonal (Ih), the bond-centered (Ib), 
and the dumbbell or split (Is) interstitials. Biswas and Hamann calculated the 
energy of formation of It, Ih, and the vacancy. Tersoff, using T2, calculate^ 
"upper bound" values for the energies of It, Ih, Ib, and the vacancy; with T3, 
he reported energies for the same defects and also for Is, lx, the split vacancy 
and the saddle point for the concerted exchange mechanism of Pandey. 
Using the procedure outlined in Sec. 3, we have performed our own calculation 
for It, Ih, Ib, Is, the vacancy, and the split vacancy. We have calculated defect 
energies at constant atomic volume (Q = fie) and at constant pressure (P - )• 
The volume relaxation reduces the formation energies by less than 0.1 eV. This 
is because our cell is very large and extends to the 24th shell from the defect. 
For a (64 ± 1) atoms cell, volume relaxation is more important and reduces the 
formation energies by up to 0.3 eV. The formation energies of the equilibrium 
(relaxed) configurations presented in Table VI include volume relaxation, i.e., 
they are formation enthalpies at zero pressure. There is no reliable experimental 
data for the equilibrium energies and structures of these defects; therefore, our 
results are only compared to those obtained with first-principles methods. As 
indicated in Sec. 3, there are, however, large uncertainties in these ab initio data 
so that only a range of values is really available at this time. These are also 

listed in Table VI. 

PTHT understimates significantly the energies of all the defects, in parti cu ar 
I T and the vacancy. The short-range potentials give a better description of the 
energy of It than the longer range potentials (PTHT and BH). This would 
confirm the observation of Biswas and Hamann that short-range functions are 
needed to model such a defect. 10 In fact the short-range potentials (T3, SW, 
DOD, and T2 in that order) appear to give a better overall description of the 

energies of all defects considered here. 

PTHT gives a small relaxation of the nearest neighbors surrounding It- With 
all potentials, It has four and six first and second neighbors at 2.44 and 2.8 A 
(PTHT), 2.54 and 2.84 A (BH), 2.56 and 2.94 A (SW), 2.57 and 2.69 A (DOD), 
2.52 and 2.72 A (T2), and 2.38 and 2.96 A (T3). Ih has six first neighbors at 
2 42 A (PTHT, DOD, and T2), 2.51 A (BH), 2.58 A (SW), and 2.48 A (T3). Ib 
has two first neighbors at 2.36 A (PTHT), 2.63 A (BH), 2.31 A (SW and DOD), 
2.27 A (T2) , and 2.23 A (T3). The two atoms forming Is are separated by 
about 2.19 A (SW) - 2.33 A (PTHT). They each have two second neighbors at 
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2.28 A (T3) - 2.36 A (SW). In general, all potentials confirm the findings of 
Batra, Abraham, and Ciraci 50 that (i) there are significant atomic relaxations 
extending to several shells around the defect and (ii) the relaxation is oscillatory 
in nature and somewhat nonuniform within some shells. 

BH and SW lead to an equilibrium configuration of the vacancy where the 

neighboring atoms relax radially (along (ill)) inward towards the defect. This 
relaxation brings the surrounding atoms (initially separated by 3.84 A) closer to 
each other to about 2.89 A resulting in a weak interaction of their dangling 
orbitals. It also increases the length of the back bonds to 2.44 A. We note that 
for SW, the ideal configuration is metastable (see Sec. 6.1.1). Biswas and 
Hamann reported that the relaxation energy for the vacancy was zero. We 
found that, unlike for SW, the ideal configuration is not at equilibrium because 
the forces on some atoms although small are quite significant. PTHT, DOD, T3, 
and T2 give an outward relaxation (also along (ill)) in analogy with the (111) 
lxl surface. In fact, the amount of relaxation correlates somewhat with the first 
interlayer contraction of that surface (see Table VIII). With PTHT, DOD, and 
T3, there is another metastable configuration for the vacancy where the 
relaxation is inward. In this case, the vacancy formation energy and the 
fractional amount of relaxation are 2.87 eV and -31.5% (PTHT), 4.21 and -32.2 
(DOD), and 4.00 and -28.9 (T3). 

Whether the relaxation around the vacancy is inward or outward is still a 
subject of controversy. Most past quantum-mechanical calculations lead to an 
outward relaxation (see Ref. 96 for a discussion of such work). More recently 
Kelly, Car, and Pantelides 96 found an inward relaxation of about 0.2 A 
compared to 0.6 A and 0.56 A for BH and SW. Antonelli and Bemholc also 
found an inward but smaller relaxation (-2.8%). Wang, Chan, and Ho using a 
tight-binding molecular-dynamics method also arrived at an inward relaxation of 
0.5 A which compares very well with the results obtained with BH and SW. 

The configuration corresponding to the split vacancy is the classical saddle 
point for vacancy migration. In this migration path, a neighboring atom moves 
along a bond and displaces the vacancy. The saddle point is expected to be the 
configuration where the atom is at the mid-bond site. This atom has six first 
neighbors at about 2.96 A when these are at their ideal positions. The ab mitio 
results 96 for this defect are: E f = 5.01 eV (unrelaxed) and 4.19 eV (relaxed), 
and an inward relaxation of the six neighbors of 0.28 A. The calculation of the 
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vacancy formation energy lead to 4.5 and 3.92 eV for the unrelaxed and relaxed 
configurations, respectively. 96 The resulting migration energy for the vacancy is 
0.27 eV compared to the experimental value of 0.45 eV. We note that both 
Tersoff potentials give the split vacancy as the most stable configuration for the 
monovacancy. The inward relaxation of the neighbouring atoms is described 
well by T3, SW, and BH. Assuming, as in Ref. 96, that the configuration of the 
split vacancy is the saddle point for vacancy migration, the migration energy for 
the vacancy as obtained with BH and SW is 0.18 and 0.54 eV, respectively. We 
note the excellent agreement with experiment provided by SW. 

Our results obtained with SW agree with those of Ref. 50. For T2, only the 
result for the vacancy agrees with that of Ref. 12. For It, Ih> and Ib, our results 
for E f are consistently smaller but, as indicated above, Tersoff stated that his 
numbers were upper bounds on Ey . Also our values for the formation energy of 
It obtained with T3 is smaller by 0.35 eV than Tersoffs result. 13 The largest 
discrepancies are with the relaxed formation energies of It and Ih obtained with 
BH. Biswas and Hamann 10 reported values of 3.61 and 5.09 eV, respectively, 
compared to our results of 1.56 and 2.89 eV. On the other hand the unrelaxed 
formation energies are in much better agreement: 4.99 and 9.47 eV compared to 
our numbers 4.57 and 9.31 eV. Our result for the unrelaxed vacancy formation 
energy is exactly the same as theirs. Biswas and Hamann did not specify clearly 
the cutoff radius and size of the cell they used for the defect calculations. To try 
to resolve this question, we performed calculations for these two defects with a 
small cell containing 65 atoms and used different values for the cutoff radius, Re- 
With R c = 5.0 A, there is very little change in the results of Table VI. With Rc = 
3.95 A, the formation energies are 0.83 eV (4.86 eV, unrelaxed) for It and 2.57 
(9.73) for Ih- Finally, with R c = 3.0 A, we found 4.35 (9.74) and 5.56 (14.46), 
respectively. We do not know at this time what is the reason behind these 
differences. We do believe that our results are correct. 

5.6 The liquid and amorphous states. 

We have not performed any simulation of liquid (/-Si) or amorphous (a-Si) 
silicon; however, numerous studies have been performed by others to study 
these two bulk phases. For a-Si, the resulting structures seem to be very much 
dependent on the simulation procedures. 52 ' 54 Also, the potentials were 
sometimes altered in order to achieve a better description of the amorphous 
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state. 43 For these reasons, we will not attempt to review the simulation work 
done with BH, SW, T2, and T3. Instead, the interested reader is directed to the 

original literature given for each potential in Sec. 2. 

Takai, Halicioglu, and Tiller performed a constant pressure Monte^Carlo 
simulation to study the melting of silicon using the PTHT potential. The 
melting temperature, T m , was determined to be about 1920 K in fair agreement 
with the experimental value of 1685 K. The potential correctly describes the 
volume contraction upon melting; however, other properties of /-Si such as the 
latent heat of fusion are 2 to 3 times smaller than the experimental values. The 

structural properties were not reported. 

Because a requirement in the fitting procedure was to accurately reproduce 

the melting point of the crystal and the structure factor of the liquid, 6 SW gives 
the best overall description of /- Si than any other potential. Several groups 
studied /-Si and a-Si using this potential and different simulation methods. 
53,55,62.67-68 stillinger and Weber determined T m to be about 2013 K. 6 However, 
all other studies lead to a value in the range 1665 - 1750 K. The structure of /-Si 
is also rather well described and other properties are in good or fair agreement 
with experiment. 

BH and T3 strongly overestimate T m which is about 2900 K ♦ and 3000 
13 , respectively. The radial distribution function of the liquid is described fairly 
well by T3. T2 does not describe /-Si well. 12 Both Tersoff potentials seem to 
favor four-fold coordination in the liquid in disagreement with experiment. 

Note that the simulations of /-Si reported here for PTHT and SW were 
performed with the original sets of energy parameters. These parameters, as 
indicated in Sec. 2, give a bulk cohesive energy, E c , of diamond silicon of -5.45 
eV for PTHT and -4.34 eV for SW. Simulations on melting have predicted strong 
correlation between T m and E c . 117 Thus we may scale T m with E c and the 
results for PTHT and SW become 1631 K and 1776 - 1867 K, respectively. 

6 SURFACES 

In view of the rich variety of surface reconstructions they exhibit, surfaces of 
silicon represent, perhaps, the most stringent test for the potentials. We consider 
in this section the low-index (100), (111), and (110) surfaces. Because of their 
technological importance in the microelectronics industry, the (100) and (111) 
surfaces of silicon have been extensively studied both experimentally and 
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theoretically over the last 30 years. A myriad of experimental techniques were 
used to determine their geometrical and electronic structure. Several models 
have been proposed and many theoretical studies, most notably self-consistent 
total-energy pseudopotential calculations, have been performed to explain the 
different surface patterns observed. For a complete review of these surfaces the 
reader is referred to the articles of Haneman. 118 It is worthwhile to note that 
experimental observations of these surfaces are made at T * 0. Temperature 
effects associated with the entropy can therefore be important and might be the 
critical factor responsible for the relative stability of some surfaces, e.g., Si(100) 
and Si(lll) 7x7. However, because calculations or even estimates of the 
entropy are difficult, we use, as is common, zero temperature surface energies 
to study the stability of such surfaces. 

6.1 Si(100) 

Despite its apparent simplicity, the Si(100) surface has long been the subject 
of controversy. Many models have been proposed over the years to explain the 
various surface patterns observed experimentally. 119 > 120 * 121 It is now universally 
accepted that the dimer model is the correct one. The dimerization of the 
Si(100) surface has been particularly confirmed by STM observations in real 
space 120 - 122 and also by numerous total-energy calculations, principally those 
using first-principles pseudopotential techniques. 86 ’ 90 * 123 126 The reduction of the 
dangling bond density is the primary driving force for the dimer reconstruction. 
Dimerization occurs when two surface atoms (initially in their ideal bulk 
positions and separated by the second neighbor distance of 3.84 A), which have 
two dangling bonds per atom, move toward each other in the [110] direction and 
in the plane containing their dangling bonds to form a bond. This is illustrated in 
Fig. 9 which shows symmetric dimers. The dimerization induces subsurface 
atomic displacements which extend at least four layers into the bulk. The dimer 
is buckled and asymmetric when the dimer atoms have different x and z 
displacements from their ideal bulk positions. The x and y directions, taken as 
[110] and [1 10], run along the dimer bond and the rows of dimers, respectively. 

Both 2x1 and c 4x2 patterns have been observed in LEED and He-atom 
diffraction experiments; the latter also showed the presence of p 2x2 and 
possibly c 2x2 patterns. 119 The situation was clarified recently by Tromp, 
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Hamers, and Demuth who used an STM with a lateral resolution of about 3 A to 
determine the atomic structure of the clean Si(100) surface. 123 They only 
observed asymmetric buckled and symmetric nonbuckled dimers along with a 
relatively high density (approximately 10%) of dimer vacancies. The surface 
defects, which are randomly distributed, appear as both individual dimer 
vacancies and small clusters of missing dimers. No other type of reconstruction 
was observed. The density of buckled and symmetric dimers is nearly the same 
indicating that their energies are approximately equal. This is supported by ab 
initio calculations as we will see below. 90 - 124 - 125 in defect-free regions, the 
dimers are symmetric and the periodicity is 2x1. The buckled dimers, which give 
rise locally to p 2x2 and c 4x2 patterns, are often observed near vacancies and 
at steps. This suggests that the surface defects induce or at least stabilize 
buckling. 120 The p 2x2 and c 4x2 structures are formed by alternating the 
buckling along a row of dimers for both of them but, in adjacent rows, the 
buckling is parallel and antiparallel for the former and latter, respectively. The c 
2x2 structure was not observed. 

In this work, the Si(100) surface is modeled with a slab containing 20 layers 
of 16 atoms each. The top 14 layers are allowed to relax while the rest are held 
fixed. The cartesian coordinate system is shown in Figs. 9 and 10. For the dimer 
reconstruction, the surface atoms are initially displaced toward each other to 
form dimers. The results for the ideal and relaxed lxl, the dimer reconstructed 
2x1 and c 2x2 surfaces, and the Pandey 7 t-bonded defect structure are 
presented and compared to the DFT results in Table VII. 

6.1.1 Si(100) lxl 

The energy of the ideal lxl surface, as obtained with the six potentials, is 
comparable to the DFT result of 2.5 eV. This value obtained at 4.3 Ry is 
probably an upper bound. The largest discrepancy is obtained with PTHT and 
DOD which give about the same value. By relaxing the lxl surface, the energy 
is slightly lowered in agreement with the DFT result. This relaxation, of both 
energy and first interlayer contraction, is perfectly described by BH and to a 
lesser extent by T3. Again, PTHT and DOD produce similar results; they both 
overestimate the relaxation energy by about a factor of 3. DOD gives an 
interlayer contraction of 10%, twice the DFT result. T2 gives very little 
relaxation. In the case of SW, all stresses are zero and the surface exhibits no 
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relaxation; this is true, as we will see later, for any bulk terminated surface or 
any defect that is created by removing atoms, e.g., vacancies. This behavior can 
be directly related to the form of the potential 37,58 and is explained by the 
combined effect of three features: (i) the potential includes only first-neighbor 
interactions, (ii) the total 3-body energy vanishes at the tetrahedral angle, and 
(iii) the Si2 molecule bond length and strength are exactly equal to the bulk 
equilibrium bond length and energy. As a result of (i) and (ii), the forces and all 
stresses on all atoms are zero; consequently, the ideal lxl structure of any 
surface is a minimum on the potential energy surface. There is no relaxation (for 
unreconstructed surfaces) because of (iii). Compared to the DFT results, the 
potentials do not, in general, predict the surface stress well. For the ideal lxl 
surface, they all underestimate o xx ; the short ranged potentials even predict a 
zero stress for that direction. PTHT and BH overestimate o yy strongly; T3 
predicts a negative value. 

6.1.2 Dimer Reconstruction 

Since we will compare our results with those obtained via ab initio methods, 
it is worthwhile to review those theoretical calculations. As stated above, 
numerous total energy calculations were performed over the years in an attempt 
to explain the experimental observations. 86 - 90 * 118 ’ 123 126 We only focus on three 
of them which all used the pseudopotential technique. 

Pandey 124 found that the lowest energy configuration in a 2x1 cell leads to 
symmetric dimers. The surface energy (relative to the ideal lxl surface), 
computed with Epw = 6 Ry, is -2.06 eV/dimer. Pandey also showed that, using 
this minimum energy structure and buckling the dimer such that no bond length 
is altered, the energy was raised by 0.02 and 0.11 eV/dimer for a tilt angle (with 
respect to the surface) of 10° and 15°, respectively. Batra 90 also obtained 
symmetric dimers for his optimized 2x1 surface with a relative energy of -2.34, - 
2.06, and -1.86 eV/dimer at 5.5, 6.5, and 7.5 Ry, respectively. Note that this 
indicates that the surface energy has not yet fully converged. An optimized 
buckled geometry with 2x 1 symmetry and a tilt angle of 8° raised the energy by 
only 0.02 eV/dimer at 7.5 Ry. However, buckled dimers in 2x2 symmetry 
lowered the energy by 0.0054 eV/dimer (this structure was not optimized and 
the buckling was very small). Batra also showed that there is no barrier to 
dimerization and that twisting of the dimers is energetically unfavorable. Roberts 
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and Needs ^ performed what is perhaps the most extensive calculation of 
dimer reconstruction on the Si(100) surface. They considered different 
geometries including the 2x1 surface with symmetric and buckled dimers, the p 
2x2 structure with alternating buckled dimers, and other structures with missing 
dimer defects. They found that buckled dimers have lower energy than 
symmetric dimers. The relative energies, obtained at 6 Ry, for the buckled p 2x2 
and 2x1 structures and the symmetric 2x1 structure are -2.108, -2.078 and -2.02 
eV /dimer, respectively. The energy difference is indeed small and is comparable 
to the accuracy of their calculation. In fact, the energy difference between the 
buckled and symmetric dimer in the 2x1 structure is only 0.03 eV/dimer at 10 
Ry. The buckled dimer in the 2x1 structure is tilted by about 6.9° while the two 
alternating buckled dimers in the 2x2 cell have different tilt angles, 11.6° and 
12.3°. For the optimized 2x1 structure, the structural parameters are nearly the 
same for the three calculations. They are given in Table VII and Fig. 9. The 
main results from these calculations are (i) the symmetric and buckled dimers 
have nearly the same energy which is compatible with the observation that they 
have nearly the same density, (ii) the dimers are not twisted, (iii) the dimer bond 
length is shorter than the bulk bond length indicating multiple bonding character, 

and (iv) the back bonds are strengthened. 

For all potentials, the Si(100) surface reconstructs to form symmetric dimers. 
Twisting the dimers raised the energy in agreement with the DFT result. No 
potential is able to model buckling which is a quantum-mechanical phenomenon. 
123,128 will probably require inclusion of atomic interactions higher than 
three-body in the potential and even possibly fitting the potential parameters to a 
buckled structure. Both the 2x1 and c 2x2 structures are found to be stable with 
a small energy difference between them. Only PTHT and DOD incorrectly 
predict that the c 2x2 surface is more stable. For PTHT, the energy gain of the c 
2x2 structure over the 2x1 structure is due to a lower two-body energy, the 
three-body energy being about the same, while for DOD, it is due to a lower 
three-body energy (positive) despite the increase of the negative two-body 
energy. The structural parameters (bond lengths and angles) are nearly the 
same for both structures, which is consistent with their small energy difference. 

We now focus on the 2x1 reconstruction. Note that, as indicated above, the 
DFT value of -0.93 eV for the relative energy of this surface is certainly a lower 
bound. Thus, all potentials, but T2, lead to a surface energy in fair agreement 
with the DFT result. Only PTHT, DOD, and T2 give a dimer bond length smaller 
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than the bulk equilibrium bond length (2.352 A). The value for T3 is very close, 
the largest discrepancy is obtained with BH and SW. For the back bonds 
between surface and second-layer atoms, the DFT calculations show a 
strengthening with a length slightly larger than the dimer bond length. Only 
PTHT, DOD, T2, and T3 give this back bond strengthening. Also, only for T2 is 
the length of these back bonds larger than that of the dimers. Note, however, 
that all bond lengths are within 2% of the bulk bond length and that the largest 
discrepancy between the DFT results and those obtained with the potentials is 
7, 3, and 5% for the dimer bond length, the length of the back bonds, and the 
bond angles, respectively. 

According to the DFT calculation, the ideal lxl surface is under a strong 
tensile stress along the eventual dimerization direction and under a weaker 
tensile stress in the perpendicular direction. The reduction of o xx after 
dimerization suggests that stress relief is the driving force for the dimer 
reconstruction. However, if this were the case, 0 yy would also be small; but it is 
now stongly compressive. Therefore, just as for the Si(lll) 7x7 structure, it is 
primarily the reduction of the dangling bond density that stabilizes the dimer 
reconstruction on the Si(100) surface. Only SW and T3 predict the correct sign 
for the two stresses. SW overestimate c xx by 68% and gives a vanishingly small 
o yy while T2 underestimate both of them by 47 and 36%, respectively. BH and 
T2 predict the value of o xx very well but both give a positive o yy (almost zero for 
BH). On the other hand, PTHT and DOD predict a yy very well; they, however, 
predict a compressive stress in the x-direction. It is worthwhile to indicate that, 
for all potentials (even those for which the sign of one of the stresses is wrong), 
the 2x1 surface is more compressed (less tensile) in the y-direction than it is in 
the x-direction. 

6.1.3 Surface Defects 

We now turn our attention to the surface defects. It does appear that the high 
defect density is intrinsic of the (100) surface. Indirect support for this is 
provided by the observation that similar sample preparation procedures for (111) 
surfaces lead to a low surface defect density 120 and by the low formation 
energy of dimer vacancies. 125 On the other hand, the geometrical arrangement 
of the defects depends on the sample preparation technique and/or the presence 
of impurities on the surface. In the STM experiments, 120 - 123 the samples were 
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gradually cooled down after annealing at high temperatures, this led to a random 
distribution of missing dimers. In another STM experiment, 129 missing dimer 
defects ordered along the dimerization direction are clearly seen when small 
amounts of gallium are deposited on the Si(100) surface. Indirect evidence for 
ordered dimer vacancies is also provided by the LEED studies of higher-order 
periodicities which can be produced by rapid quenching from high temperatures. 
These are c 4x4 and c 8x8 130 and, in particular, 2xn with 6 < n < 10. 131 These 
structures have all been explained with models involving ordered dimer 
vacancies. For the 2xn structures, it has been determined that they are 
metastable and that quenching at higher rates or from higher temperatures lead 
to small n. 131 It is clear from the above that dimer vacancies play an important 
role in the energetics and structures of the Si(100) surface. For the potentials to 
be useful, it is essential that they give at least a reasonable description of such 
defects. Roberts and Needs 125 calculated the surface energy of the Pandey n- 
bonded defect structure. 124 In this model, shown in Fig. 10, a dimer vacancy is 
created at every fourth site along the rows of dimers. This leads to a structure 
with 2x4 symmetry. In the DFT calculation, the structure was optimized but the 
dimers were not allowed to buckle. Its energy is only 0.035 eV/lxl cell higher 
than that of the 2x1 symmetric dimer reconstruction. This corresponds to a 
formation energy (with respect to the 2x1 symmetric dimer structure) of 0.28 
eV/defect. In a 2x2 cell, where the dimer vacancy concentration is 0.5, the 
formation energy is 1.10 eV/defect. Clearly these surface defects repel each 
other. Now, while the creation of a dimer vacancy induces strain in the surface, 
it does reduce the dangling bond density by rebonding of the exposed second- 
layer atoms. Thus, the optimum surface defect concentration is probably lower 
than 0.25 and is a compromise between the two effects. 

The results of our own calculation for this defect structure are presented in 
Table VII. DOD considerably overestimates the strain energy and leads to a 
surface energy only slightly lower than the ideal lxl surface. At the other end, 
T2 predicts that vacancy formation is exothermic (with respect to the 2x1 
structure) in agreement with Pandey ’s predictions (Pandey 124 , who did not 
perform a calculation, estimated that the defect would lead to a large energy 
gain of 2.0 eV/defect which seems quite unreasonable). BH , SW, and T3 do 
predict that the defect structure is metastable and give results in fair agreement 
with those of the DFT calculations. While, as stated above, the potentials do not, 
in general, give a good description of the surface stress, it is worthwhile noting 
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that, qualitatively, the surface defects induce atomic displacements such that the 
tensile stress along the dimer bonds is reduced and that the surface is now 
under tension (or at least under less compression) in the y-direction. The atomic 
displacements are qualitatively similar to those of the DFT optimized structure. 
Only DOD predicts a shortening of the bond length of the dimers adjacent to the 
defect; the others predict a small expansion (< 1%). The length of the weak 
bond formed by the exposed second-layer atoms is 2.562, 2.625, 2.318, 2.871, 
and 2.656 A for BH, SW, DOD, T2, and T3, respectively compared to the DFT 

result of 2.71 A. 

In conclusion, BH, SW, T3, and to a lesser extent T2 give a reasonable 
description of the energetics and structures of the Si(100) surface whose 
principal features are the dimerization of the surface atoms and the existence of 
dimer vacancies. This is mainly due to the fact that the angular distortions (from 
the tetrahedral angle) on this surface are relatively small (about ± 12% for the 
2x1 surface; they are somewhat larger for the defect structures), at least 
compared to those encountered in small clusters and on the (111) surface. Their 
principal limitation is the inability to model buckling which might not be serious 
because the energies of the buckled and symmetric dimers are nearly 
degenerate. Clearly surface defects play an important role here and any realistic 
calculation or simulation which aims to study them will have to involve systems 
with a large number of atoms. Powerful and accurate first-principles methods, 
such as the now widely-used ab initio molecular dynamics technique, are still 
limited to studies of relatively small systems (up to 150 atoms 133 ). It is believed 
that the potentials mentioned above will be useful in large-scale simulations of 
the Si(100) surface since they predict its properties reasonably well and they 
can handle large systems with atoms in the tens of thousands. 60 

6.2 Si(lll) 

The Si(lll) surface exhibits several reconstructions depending on the purity 
of the surface, the temperature, and the sample preparation procedures. We 
consider here only clean surfaces. A freshly cleaved (111) surface in UHV 
reconstructs to form a 2x1 metastable structure which transforms irreversibly 
upon heating (the transition temperature ranges from 200 to 350 °C depending 
on the step density 134 ) to the stable 7x7 phase which, in turn, transforms 
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reversibly to a structure with lxl symmetry at about 830 °C. This lxl phase, 
which still remains a mystery, 118 is not the bulk terminated surface. 

The widely accepted model for the 2x1 surface is the Pandey u bonde 
chain model. 85 > 91 ’ 132 ’ 135 ' 138 This structure is similar to the (110) surface. Another 
interesting model, called the three-bond-scission model, has been recently 
proposed by Haneman. 118 ’ 139 For the 7x7 surface the DAS model of 
Takayanagi, Tanishiro, Takahashi, and Takahashi is now universally accepted. 
mo- 144 other metastable structures have also been observed. STM expenments 
on (111) surfaces prepared by a combination of laser and heat treatment 
showed locally the existence of 2x2, c 4x2, ^3x^3, 5x5, and 9x9 structures. 

The first three have been explained by simple adatom covering of the surface 
with the adatom located in the T 4 site. The (2n+l)x(2n+l) structures were 
explained with DAS type models. 144 ' 146 A ^3 structure has also been 
observed in a LEED experiment; 147 the LEED data were fitted to a vacancy 
model which was shown to be energetically unfavorable. 

We have performed calculations for most of the structures mentioned above. 
In general, the Si(lll) surface was modeled with a slab containing eight double 
layers. The top four or five double layers were allowed to relax while the rest 

were held fixed. The orthogonal x and y surface axes run along the [110] and 
[1 12] directions, respectively. Depending on the symmetry of the surface, the 
number of atoms per layer ranged from 24 to 162 corresponding to a total 
number of moving atoms of 240 to 1640, respectively. The surface energies and 
stresses for the (111) surfaces are presented in Table VIII. 

6.2.1 Si(lll) lxl 

Compared to the DFT result, the surface energy of the ideal lxl surface is 
underestimated by all potentials. Although bond breaking energies are certainly 
important, only the relative energies are relevant when investigating the stability 
of one surface structure over another. We note that the relaxation of the lxl 
surface is best described by PTHT, DOD and T3; however, they all significantly 
overestimate the compressive lateral surface stress. With BH and T2, the lxl 
surface exhibits little relaxation resulting in a small tensile stress. As indicated 
above, in the case of SW, all stresses are zero and the surface exhibits no 
relaxation. As expected, all short-ranged potentials predict a zero surface fault 
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energy for the lxl faulted surface. PTHT and BH give a value of 0.012 and 
0.015 eV/lxl cell compared to the DFT result of 0.02 to 0.06 eV. 91 They, 
however, predict a surface stress of -2.31 and -0.04 eV/lxl cell while the DFT 
value is 0.11 eV/lxl cell. 

6.2.2 Si(ll 1)2x1 

Using a cleavage technique, Gilman ^ measured the surface energy of the 
(111) surface at 77 K. He obtained a value of about 1 eV/lxl cell. Since, as 
indicated above, a freshly cleaved surface transforms to a 2x1 reconstruction, 
except for any generated strain energy associated with the cleavage process, 
this value can thus be taken as the experimental surface energy of the (111) 2x1 
7 t-bonded surface. The DFT result of 1.11 eV compares very well with this 

value. 

Only PTHT predicts that the 2x1 7 t-bonded structure is stable with respect to 
the ideal lxl surface, the relative energy is however underestimated by a factor 
of 3. Moreover, the DFT calculations predict that the structure is also buckled. 
132,135,137,138 Again the potentials cannot model buckling. For BH, this structure 
is not even a minimum. The 2x1 surface reduces its energy by 7 t-bonding of the 
surface dangling bonds which are now first neighbors instead of second 
neighbors as they are on the lxl surface. Since the potentials do not model 7t 
bonding, for them the dangling bond density is unchanged. Therefore, there is no 
energy gain to overcome the strain energy caused mainly by angular distortions 
(±10%). Only for PTHT is this strain energy compensated by a strengthening of 
the bonds in the surface chains. DOD shows a somewhat similar behavior, but 
the strain energy is higher due to a stronger three-body energy (Fig. 4). Our 
result for DOD is in disagreement with that of Dodson who reported that the 2x1 
reconstruction reduces the energy of the (111) surface by 0.12 eV/lxl cell. In 
general, the potentials predict that the 2x1 surface is under a weak and stronger 
tensile stress in the directions parallel and perpendicular to the chains , 
respectively; the DFT calculation predicts the opposite. 85 However, the 
comparison here is not truely appropriate because the DFT structure is buckled 
and stresses are more sensitive to the atomic displacements than is the energy. 

6.2.3 Adatom Structures 
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According to the DFT calculations, 91 all adatom covered structures are 
stable (with respect to the ideal lxl surface) in the order 2x2 T4 < ~'l3x y l3 T4 < 
2x2 H3 < V3x^3 H3. Thus, the 2x2 structure is more stable than the V3xV3 
structure in spite of the fact that it has a smaller reduction of the dangling bond 
density. The energy difference is, however, only 0.06 eV/lxl cell. Meade and 
Vanderbilt have attributed the energy lowering to the different electronic 
structure of the two surfaces. 91 Note that, with the same adatom concentration 
of the 2x2 structure, two different structures with rectangular c 2x4 and c 2x8 
symmetries can be generated. 146 We have not considered these structures in 
this work. Before we discuss the results for the adatom structures, note that 
there are four errors in Table I of Ref. 37. For T3, the relative energy of the 
relaxed lxl surface should be -0.07 and the energy and stress of the V3xV3 H3 
surface must be 0.482 and -0.522, respectively. The surface stress of this latter 
structure should read -0.502 for BH. Also, our calculation with SW for the 
y l3x y l3 structures do not agree with that of Ref. 27. Li, Chen, Allen, and 
Broughton reported 1.53 and 1.82 eV/lxl cell for the surface energy of the 
V3xV3 H3 and T4 structures; 27 our numbers are 1.12 and 1.61. 

Only T2 predicts that all adatom structures are stable as they should be; 
however, the relative energies are strongly underestimated. With BH, DOD, and 
T3 all adatom structures are unstable. PTHT predicts that only the H3 
structures are stable. For SW, only the V3xV3 H3 surface is stable. With the 
exception of DOD, all potentials incorrectly favor close packing of adatoms. 

According to the DFT calculations, 91,93 in both the T4 and H3 structures the 
three upper atoms surrounding the adatoms relax laterally inward. In the T4 
structures, the adatom is bonded to the three upper atoms and to the second- 
layer atom directly below; the bond lengths are about the same for both types of 
neighbors. It is 2.47 and 2.49 A in the 2x2 91 and ^3xV3 93 structures, 
respectively; the bond angles are about 92® and 56®. In the H3 structures, only 
the three upper atoms are first neighbors of the adatom; for the V3xV3 structure 
the bond length is 2.55 A and the bond angle is 93®. 93 The inward relaxation of 
the surface atoms in the T4 structures is predicted by all potentials while in the 
H 3 structures, only BH, SW, T2, and T3 predict such relaxation. BH and, in 
particular, SW overestimate significantly the separation between the adatoms 
and their neighbors in the T 4 structures. For T2, the atomic displacements in the 
structures are in reasonable agreement with those obtained in the DFT 
calculations. In general, for all potentials, the structural parameters of the T4 
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structures are qualitatively in better agreement with the DFT results than those 
of the H 3 structures. Consequently, because the surface stress is strongly 
influenced by the atomic displacements it is qualitatively better described in die 
T 4 structures. These adatom structures are under a strong tensile stress which 
is consistent with the inward relaxation of the three surface atoms surrounding 
the adatom. The greater this relaxation, the larger the positive surface stress. 

For all potentials with the exception of T2, the energy gain resulting from the 
reduction of the dangling bond density is not enough to overcome the strain 
energy caused primarily by the very small bond angles (Fig. 4). For T2, because 
the angular function is more flexible and the bond bending forces are small (see 
Sec. 2), it is able to give a better description of the energetics of these surfaces 
than the other potentials. Even for those cases where the relative energy of the 
various adatom structures is positive, the adsorption energy is negative as it 

should be. 

6.2.4 Si(lll) V3x^3 - Vacancy Model 

For the vacancy model of the ^3x^3 surface, all potentials correctly predict 
that this structure is not stable. The vacancy formation energy (with respect to 
the ideal lxl surface) is however significantly overestimated by all potentials 
with the exception of SW which predicts a value of 0.24 eV/vacancy compared 
to the DFT value of 0.42 eV/vacancy. The DFT calculation predicts a moderate 
inward relaxation (toward the vacancy) of the second-layer atoms and a large 
compression of the first interlayer spacing. No value for the surface stress is 
available. With BH and SW, the structure is under a strong tensile stress 
reflecting the large inward relaxation of the second layer atoms which results in 
a large bond stretching. On the other hand, PTHT, DOD, and T3 give a 
compressive stress which also reflects the outward relaxation in this case. 
Finally, consistent with the lack of relaxation, T2 gives a very small stress. For a 
complete discussion of the vacancy structure the reader is directed to Refs. 37 

and 148. 

6.2.5 Si(lll) (2n+l)x(2n+l) 

The DAS model for the Si(lll) (2n+l)x(2n+l) surface is well known and the 
reader is directed to Refs. 141, 144, and 145 for a complete description of the 
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structure. This surface contains several structural units which are the dimers 
along the domain walls, T 4 adatoms with a local 2x2 symmetry, stacking faults, 
and corner holes (an extended surface vacancy). The surface energy is the 
result of a balance between the individual contributions of these units and 
possibly the interactions between them. This surface is obviously the most 
stringent test for the potentials. 

Our results, presented in Table VIII, show that only T2 predicts that the 7x7 
DAS structure is stable with respect to the ideal lxl surface. Furthermore, it is 
the ground state of the ( 111 ) surface at least compared with the structures 
considered so far. However, as before, the energy is underestimated. 
Considering the fact that it is also only T2 which predicts that the 2x2 T 4 
structure is stable, we conclude that the adatoms play a major role in 
determining the energy of the 7x7 surface. To confirm this, we have also 
calculated the energy of the same surface but without the adatoms, the so- 
called DS model . 85 The relative energy of this structure is indeed very small; it 
is negative in the case of PTHT, DOD, and T2. The energy difference between 
the DAS and DS structures should give the contribution of the adatoms to the 
total energy of the DAS structure if their interaction with the other structural 
units is negligible. This energy difference is 1.47, 1.25, 1.92, 1.69, -0.48, and 2.39 
eV/adatom for PTHT, BH, SW, DOD, T2, and T3, respectively. Except for 
PTHT, this energy compares fairly well with the contribution of the adatoms to 
the 2x2 T 4 structure (0.93, 1.38, 2.05, 1.82, -0.31, 2.54 eV/adatom) indicating 
perhaps that the adatoms do not interact (or very little) with the other features 
on the surface. Note that the energies of the faulted and unfaulted 2x2 T4 
structures are nearly identical. 85 91 In the representation provided by these 
potentials, the adatoms play an important role in determining the energy of the 
7x7 surface. All potentials predict that the 7x7 DAS surface is under tension in 
agreement with the estimate of Vanderbilt . 85 Considering the small value of the 
stress of the corresponding DS structure, most of the tension is caused by the 
adatoms. Our results obtained with SW and T2 for the 7x7 DAS surface agree 
with those of Ref. 27 where the structural parameters of the optimized 
structures were also reported. The investigation of the vibrational spectrum of 
this surface showed that, despite their limitations, these two potentials are able 
to accurately describe the z-polarized adatom vibrations . 22 

It has been shown that the different (2n+l)x(2n+l) structures are very close 
in energy. 16 ’ 85 ’ 140 For example, Qian and Chadi, using a tight-binding method, 
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found that the 7x7 DAS structure is only 0.008 eV/lxl cell lower in energy than 
the 5x5 DAS structure. 140 It is natural to question whether T2 is able to predict 
such a trend and also whether the 7x7 reconstruction is indeed the lowest- 
energy structure. Thus, we have performed calculations for (2n+l)x(2n+l) DAS 
and DS structures with n = 1 - 4. The relative energy and surface stress for 
these structures are plotted as a function of n in Fig. 11. T2 does predict the low 
energy difference between the different structures, about 0.1 eV/ lxl cell; 
however, contrary to what was previuosly thought, 27 the 3x3 DAS structure is 
the ground state of the Si(lll) surface in contradiction with the experimental 
observation. The energies of the DS and DAS structures appear to increase 
linearly with n for n > 2. The surface stress, on the other hand, decreases with 
increasing n. 

We can estimate the energies of the different surface structural units in the 
DS structures by using Vanderbilt’s model based on non-interacting surface 
units. 16 - 144 ’ 146 In this model, the surface energy per lxl cell of the 
(2n+l)x(2n+l) DS structures, Ay^s* relative to that of the relaxed lxl surface, 

p, is given by, 16 


^Yds 


n(2n + l)Af + 2nAw + Ac 
(2n + l) 2 


( 8 ) 


3 

where Af = y^f 1 , . -p is the lxl surface faulting energy; Aw = -d-p is the 
'lxl -faulted v l 

relative domain wall creation energy (d is the dimer energy); Ac = c — p is the 

relative corner hole energy. Using p = 0.702 eV/lxl cell, a least square fit yields 

Aw = -0.45 eV (d = 0.17 eV), Ac = 0.37 eV (c = 1.07 eV), and Af = 0.0009 eV. 

The fact that the model gives a vanishingly small value for Af (recall that T2 

predicts that Af = 0) proves that the assumption of non-interacting structural 

units is essentially valid. In Ref. 146, the parameters used in (8) (obtained from a 

combination of DFT and Keating potential calculations) are Aw = -0.66 eV (d = 

0.53 eV), Ac = 1.26 eV (c = 2.71 eV), and Af = 0.06 eV. T2 appears to give too 

small a value for the relative corner hole energy. 

Similar calculations using the other five potentials show that the surface 
energies and stresses of the (2n+l)x(2n+l) DS and DAS structures decrease 
with increasing n. For n = 4, the surface energies are already almost equal to 
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those of the relaxed lxl and the 2x2 T 4 structures, respectively. The energies of 
the DS structures were also fitted to (8) using p and Af as given in (6.2.1) and 
Table VIII. When Af was left as a free parameter, the resulting faulting energy 
was in poor agreement with the actual value. The results obtained with PTHT, 
BH, SW, DOD, and T3 are: Aw = -0.32, 0.68, 0.34, 0.62, and 0.68 eV and Ac = 
3.15, 0.82, 0.97, 1.63, and 1.27 eV. Like T2, PTHT gives a negative relative 
domain wall creation energy. The other four potentials give a relative comer 
hole energy that is in fair agreement with the result of Ref. 146. These potentials 
favor large n because Ac is too high (PTHT) or Aw is positive (BH, SW, DOD, 
and T3). 

6.3 Si(110) 

Like the (100) and (111) surfaces, the Si(110) surface exhibits various 
reconstructions; 151 however, in contrast to the former, very little interest has 
been paid to this surface. Only the ideal and relaxed lxl surfaces are 
considered here. No attempt was made to look for reconstructed structures. 
The bulk terminated surface is formed by chains of atoms running parallel to 

[110] (chosen as the x direction; the y direction is taken as [00l]). Like the (111) 
surface, this surface has one dangling bond per surface atom pointing in the 

(111) direction. Adjacent atoms along a surface chain have dangling bond 

alternating in directions (pointing in the [1 1 1] and [1 1 1] directions) and making an 
angle of ± 35.3° with the surface normal (they are all normal to the surface on 
the (111) surface). This surface is modeled with a slab containing 15 layers of 
24 atoms each. The bottom four layers were held fixed. The surface energy and 
stresses are presented in Table IX. 

Note that, because they have the same number of dangling bonds per 
surface atom, the ideal (111) and (110) lxl surfaces have also the same surface 
energy when it is expressed in eV/lxl cell. For the long-range potentials, PTHT 
and BH, y is slightly different because atoms in the top layers have slightly 
different number of neighbors on the two surfaces. The same is true for the 
surface stress along [lTO] . Recall that the surface stress of the (111) lxl surface 
is isotropic. As for the (100) and (111) surfaces, PTHT, DOD, and T3 give a 
larger relaxation energy (which correlates more or less with the first interlayer 
contraction) than BH and T2. As discussed earlier, there is no relaxation with 
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SW. Unlike the (100) and (111) surfaces where relaxation involved only normal 
displacements of atoms (mainly of the top two layers), the surface atoms of the 
(110) surface (and to lesser extent those of the second layer), relax radially 

inward along the (111) direction. That is, in addition to the inward displacement 
normal to the surface, there is a smaller displacement in the plane of the surface 
and perpendicular to the chains. This lateral displacement results in a shortening 
and an increase of the bond lengths and angles in the chains, respectively (2-29 
A and 114° for PTHT and T3, 2.31 A and 112° for BH and T2, and 2.26 A and 
116° for DOD). The trend in surface energy for the lxl surfaces is Y(lll) < 
Y( 1 1 0) < Y(100) with all potentials. The same trend is obtained for the relaxation 
energy, Ay, with PTHT, DOD, and T3. With BH and T2 it is Ay( 100) < A Y(110) < 
AY(ill)and Ay(l 10) < Ay( 11 1) < AY(100), respectively. Note that for this 
comparison these energies were first converted to the true units of surface 
energy, e.g., eV/A 2 . 


7 OTHER POTENTIALS 

Recently, several new potentials have been proposed. These models, either 
inspired by earlier attempts or using new schemes, were intended to overcome 
the limitations of their precursors, i.e., the potentials considered in this work. It is 
thus worthwhile to review some of them. We should also mention a new class of 
total-energy functionals for semiconductors which are based on an approximate 
quantum-mechanical analysis. 32,152 

7.1 Kaxiras and Pandey 

Kaxiras and Pandey 21 constructed a potential, very similar in form to BH, in 
order to specifically simulate processes in the bulk diamond lattice. The potential 
was fitted to the entire energy surface of atomic exchange obtained from an 
accurate DFT calculation. 115 It correctly predicts the static properties of the 
perfect diamond lattice and reproduces the energy of the concerted exchange 
path to better than 0.1 eV. However, the energies of bulk point defects in their 
unrelaxed configuration appear to be too low. For the high-coordination crystal 
phases, the results are qualitatively similar to those obtained with BH. The 
potential was not tested for surfaces and clusters but it is expected that its 
predictions would also be similar to those of BH. Because the potential describes 


48 


a large range of local distortions from the perfect tetrahedral configuration very 
well, it should be useful in simulations of systems such as amorphous structures 
where the coordination remains predominantly four-fold. 

7.2 Mistriotis, Flytzanis, and Farantos 

Mistriotis, Flytzanis, and Farantos 19 modified the SW potential in order to 
correctly describe clusters with more than six atoms. The angular dependence 
of the three-body term was modified and they added a four-body term. The 
modified potential has not been extensively characterized. It has not been tested 
for surfaces. It predicts Tm to be about 2050 K and the high-density crystal 
phases are not well described. 

7.3 Khor and Das Sarma 

Following Abell 153 and then Tersoff, 8,12 Khor and Das Sarma 14 developed a 
universal interatomic potential for tetrahedrally bonded semiconductors. The 
original potential for silicon gives an excellent description of the static properties 
of cubic diamond as well as the other high-density crystal phases. The potential 
had to be somewhat extended to correctly treat surfaces. 15 The bond-bending 
term was modified to deal with the larger angular distortions from the 
tetrahedral angle encountered on the various (111) surfaces. Also, because the 
bonds of a given atom can be of a different nature, they had to fix, in an ad-hoc 
manner, the value of the effective coordination number and assumed that the 
character of the bonds remains unchanged in the course of the simulation. The 
modified potential gave a good description of the various (111) surfaces 
including the adatom structures and the Pandey 7 t-bonded chain model for the 
2x1 surface. The results for the (100) surface are similar to those predicted by 
T2. Yet another modification of the bond-bending term had to be made in order 
to successfully model the (111) (2n+l)x(2n+l) DAS surfaces. 16 It is not clear 
how these modifications affect properties previously determined. The original 
potential and its two modifications are in fact three different potentials just as 
are T2 and T3. Finally, these potentials have not been tested for bulk point 
defects and small clusters. 

7.4 Bolding and Andersen. 
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Bolding and Andersen 20 developed a potential which is a generalization of 
the Tersoff potential. The attractive term is expressed as a sum of c- and n- 
bonding terms. Interactions up to five-body are included in the potential. The 
functional form is complicated and there are over 30 parameters in this potential. 
The fitting data base was very large and included the static properties of the 
cubic diamond phase, the fact that the first pressure-induced phase 
transformation from cubic diamond is to the f3-tin phase, and finally the energies 
and geometries of global and local minima for clusters of 2 - 10 atoms. For small 
clusters, this potential generates a surface that has most of the minima (global 
and local) predicted by the ab initio calculations. However, there is not, in 
general, a one-to-one correspondence between the minima. The ground state 
structures are predicted to have energies that are in excellent agreement with 
those of quantum-mechanical calculations. The static properties of cubic 
diamond silicon are well described but the potential fails to predict the negative 
Cauchy discrepancy. For bulk point defects, only the vacancy is predicted to 
have a formation energy in good agreement with the DFT results. The energies 
of the interstitials, in particular the tetrahedral interstitial, are underestimated. 
The (111) 2x1 surface is well described. For the adatom structures of the (111) 
surface, this potential predicts results that are qualitatively similar to those 
obtained with T2. However, unlike T2, it predicts that the (111) 7x7 DAS 
surface is unstable with respect to the ideal (111) surface. Finally, for the (100) 
surface, the predictions are similar to those of T3. 

7.5 Thermodynamic Interatomic Force Field 

Chelikowsky, Phillips, Kamal, and Strauss 18 developed an interatomic 
potential similar in form to Tersoff’s. The motivation for constructing this 
potential was to study the metallic to covalent transition which occurs in clusters 
when the cluster size reaches a critical size. The angular dependence of the 
bond-bending forces was intended to describe such a transition. The potential 
describes the perfect diamond structure and the high-density polymorphs of 
silicon very well. To model clusters, it was found necessary to introduce a so- 
called dangling-bond vector which describes the transfer of bond strength from 
a dangling bond to back bonds. The energies of Si n clusters with n<10 are, 
however, still underestimated. Also, the ground state structures are in general 
not correct , e.g., for n = 3, 4 and 6. The authors predict that their potential 
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should be more useful for n > 10. The potential was not tested for bulk point 
defects and for surfaces. 

7.6 Modified Embedded Atom Method 

Baskes, Nelson, and Wright 11 proposed a new potential based on the 
embedded atom method. 4 The modification consisted of the introduction of an 
angular dependence in the host electron density. This was necessary for an 
adequate description of the bond- bending forces in the diamond cubic structure. 
This potential was extensively tested by its authors. It gives a fit to the energies 
of the high-density polymorphs comparable to some of the potentials considered 
here. It describes exactly the static properties of cubic diamond. In particular, 
unlike most other potentials it does predict the negative Cauchy discrepency. It 
also provides a fair description of bulk point defects (in particular the vacancy); 
however, it gives a high value for the energy of the intrinsic stacking fault in 
silicon. In general, surfaces are poorly described. The potential predicts an 
outward relaxation of the surface layer for all surfaces. The description of small 
clusters is in general poor. 

8 DISCUSSIONS AND CONCLUSIONS 

We have performed extensive calculations on clusters, bulk phases, and 
surfaces using the PTHT, BH, SW, DOD, T2, and T3 potentials for silicon. In 
general, no potential is able to model properly all the equilibrium structures and 
energies of small Sin (n = 2 - 6) clusters. More importantly, they all predict 
many spurious minima on the potential energy surface of these microclusters. 
The potentials do, however, predict, like the ab initio calculations, that the 
structures derived from crystal fragments are not energetically favorable even 
though the potentials were built from crystal data. T2 gives the best overall 
description of these small clusters. 

In general with the exception of the Dodson potential, the potentials do not 
accurately describe the energies of the high-pressure bulk phases. A two- 
dimensional structure with hexagonal symmetry is predicted by PTHT as the 
most stable structure instead of the diamond cubic phase. Only T3 correctly 
predicts the first pressure-induced phase transformation from diamond cubic to 
the (5-tin phase with transition pressure and volumes which are in excellent 
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agreement with experiment. SW, T3, and to a lesser extent DOD, describe the 
elastic properties well. T2 also does a good job with the exception of the 
vanishingly small value of C44. BH gives the best description of the phonon 
frequencies even though it overestimates the elastic constants. 

T3, SW, DOD, T2, and BH, in that order, give a fair overall description of the 
structures and energetics of intrinsic defects. They should be useful in studies of 
extended defects. PTHT underestimates strongly the energies of these defects. 
T2 and BH also underestimates significantly the energies of the split vacancy 
and of the tetrahedral interstitial, respectively. Only BH and SW correctly 
predict the apparently now more accepted inward relaxation of the neighbouring 
atoms surrounding the vacancy. SW yields a migration energy for the vacancy 
which is in excellent agreement with experiment. 

BH, SW, T3, and to a lesser extent T2 should also be useful in large-scale 
simulations involving the (100) surface because their predictions of its energetics 
and structures are in good agreement with those of the first principles 
calculations. On the other hand, none of the potentials is able to model the 
various reconstructions of the (111) surface. None of them model the 2x1 
reconstruction of the (111) surface correctly because of their inability to model 
7t-bonding which stabilizes this reconstruction. PTHT predicts that this 
reconstruction is stable with respect to the ideal lxl surface but with a relative 
energy that is too small. T2 predicts that the ground state for the (111) surface 
is the 3x3 DAS structure instead of the 7x7 DAS structure as previously 
thought. As for the microclusters, T2 gives perhaps the best overall description 
of the (111) surface. 

BH and SW tend to overestimate bond lengths. In fact these two potentials, 
along with T3, have, in general, similar predictions, and so do the PTHT and 
Dodson potentials. These similar behaviors correlate with similar angular 
variations of the three-body potentials, in form not necessarily in strength. The 
difference in strength is somewhat compensated for by an opposite difference in 
strength of the two-body potentials. T2 is dissimilar precisely because its three- 
body potential differs markedly from that of the others. That these similarities 
exist is quite remarkable since, except for SW and BH, these potentials are quite 
different in schemes, functional forms, and range of interactions. This attests to 
the importance of the bond-bending forces in these low-order potentials. Besides 
the fact that they do not model 71-bonding, the main reason behind their inability 
to be more transferable is an inadequate description of the angular forces. They 
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either favor only small angular distortions around the tetrahedral angle (BH and 
SW, around 127° for T3) or configurations where bond angles are much larger 
than the tetrahedral angle (PTHT and DOD). They all penalize, with the 
exception of T2, angles smaller than about 90° when in fact many structures in 
silicon involve such small angles, e.g., microclusters and some (111) surfaces^ 
The angular function of T2 is more flexible. It appears that each type of 
environment, i.e., bulk, clusters, and perhaps surfaces, needs its own angular 
function. The combined function should then be oscillatory in nature and would 
be determined by an appropriate selection of relevant structures and energies in 
the fitting database. Bolding and Andersen did just that. 20 The resulting potential 
models clusters rather well and gives a description of bulk properties 
comparable to that of some of the potentials considered here. However, it still 
failed to model properly (111) surfaces with the exception of the 2x1 
reconstruction despite the fact interactions up to five-body were included in the 
potential. This leads us to believe that it is perhaps not possible to construct a 
totally global or transferrable potential. 

In conclusion, none of the potentials considered in this work appear to be 
superior to the others. Each has its strengths and limitations. None is totally 
transferrable. Despite their shortcomings, we do believe that, some of these 
potentials will be useful in large scale simulations of materials-related problems 
as they can give insights into phenomena which are otherwise intractable to 
investigate either experimentally or with first principles methods. 
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TABLE I. Parameters for the potentials. The units are such that energy and length 
are in eV and A, respectively. Go is in degrees. 



PTHT 

BH 

sw 

DOD 

T2 

T3 

Rc 

7.3 

5.0 

3.77118 

3.2 

3.2 

3.0 

M- 


0.312058 

2.0951 

0.4 

0.4 

0.3 

G 


3.9527357 





Ai 

50928.2584 

142.2922 

189.360881 

1614.6 

3264.7 

1830.8 

A 2 

697.005028 

107.0338 

16.31972277 

155.08 

95.373 

471.18 


12 

0.5200836 

4 

2.7793 

3.2394 

2.4799 

X2 

6 

0.4206931 

0 

1.3969 

1.3258 

1.7322 

Z or Z] 

2949.47219 

26.0598 

48.61499998 




Z 2 


1.3441478 





a or ai 


0.3034373 

1.2 

4 

1.3258 

1.7322 

OC2 


0.3191903 





n 




0.6207 

22.956 

0.78734 

P 




0.13420598 

0.33675 

1.0999xl0' 6 

n 




0.8543 

2.80755739 

105.2851343 

8 




3.9588 

2.0417 

16.218 

Go 


109.471221 

109.471221 


90 

126.745381 


64 


TABLE II. Equilibrium properties of Sb and Si3. r e and r (A) are the bond lengths, 
D e and E B (eV) are the binding energies, co 0 (cm- 1 ) is the vibrational frequency, 
and 6 (degrees) is the bond angle. Values in parenthesis correspond to a 
mechanically stable configuration lying higher in energy. In the second column, the 
experimental values are for Si 2 (Ref. 97) and KR corresponds to the ab initio results 

for Si 3 (Ref. 88). 


Experiment/KR 

PTHT 

BH 

D e 

3.24 

2.38 

2.49 

re 

2.246 

2.295 

2.233 

CD 0 

510.98 

794 

463 

Eb 

7.7 
(~ 7.6) 

5.29 

(5.23) 

5.39 

e 

77.8 

(60) 

180 

(60) 

60 

r 

2.17 

(2.26) 

2.27 

(2.39) 

2.42 


SW 

DOD 

T2 

T3 

Si2 

2.32 

3.61 

2.62 

2.67 

2.352 

2.192 

2.313 

2.295 

462 

521 

467 

471 

Si 3 

4.74 

(4.63) 

7.04 

(6.05) 

7.87 

5.33 

(4.81) 

60 

(109.47) 

180 

(60) 

60 

126.75 

(60) 

2.56 

(2.35) 

2.20 

(2.40) 

2.31 

2.30 

(2.50) 
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tart f ITT Pronerties of Su - Sk. The first three most stable structures are given in the order of 
decreasing binding energy. For each structure, the first entry is the binding energy Eg (eV), others are 
tend knlthTS) and angle (degrees). KR corresponds to the ab initio results of Ref. 88. The asterisk 
indicates fhat this structure is known to be not a minima on the quanwm-m^hamcal surface The i e 

structures along with the corresponding structural parameters are illustrated m Fig. .5. lne acm^ms • 
t^Son (T d y edge-capped rhombus (ECR); corner-capped mangle and rhombus (CCT.CCR), 
pemagonal, squari, and dis.of.ed pentagonal pyramid (PP.SP.DPP); flat ^ ttogaKd £^^amrd 
(FTB ETB); orthorhombic bipyramid (OB); edge-and face-capped tngonal b.pyramrd (ECTBf CTB). For 
the linear structures, the numbering of atoms is from one end to the other. 


KR 


PTHT 


BH 


SW 


DOD 


T2 


T3 


Si4 rhombus 

rhombus 

square* 

12.85 

8.48 

8.93 

ri2 = 2.30 

2.37 

2.34 

ri3 = 2.40 

2.59 

Td 

D2d 

square* 

11.53 

8.44 

8.25 


2.34 

2.56 

Td 

CCT* 

CCT* 

9.71 

8.33 

7.38 

2.46 

ri2 = 2.26 

r 12 = 2.30 


ri3 = 2.36 

ri3 = 2.45 

Sis FTB 

ECR 

pentagon 

16.70 

12.04 

12.08 

r ]2 = 3.26 

ri2 = 2.34 

2.28 

ri4 = 2.34 

r]4 = 2.56 


r 45 = 2.78 

r 24 = 2.37 


ETB 

pentagon 

ESP* 

15.62 

11.88 

11.90 

2.48 

2.31 

2.47 

2.40 

3.86 


2.64 

FSP 

CCR 

FTB 

13.90 

11.65 

11.82 

ri 2 = 2.48 

ri2 = 2.34 

3.02 

ri5 = 2.36 

r 2 4 = 2.59 

2.46 

ri5 = 2.26 

3.48 

Si6 ECTB 

hexagon 

wedge 

21.91 

15.31 

15.81 


2.28 

n 2 = 2.53 
r]4 = 2.40 

hexagonal chair 

PP 

PP 

15.98 

14.72 

15.71 


ri2 = 2.41 

2.41 


ri6 = 2.50 

2.63 


asymmetric 

asymmetric 


14.63 

15.46 


ri2 = 2.47 

3.26 


r34 = 2.63 

3.63 


square* 

8.69 

2.39 

Td 

7.13 

2.72 


chain 

6.95 

2.35 

109.5 

pentagon 

11.57 

2.35 


FTB 

11.46 

3.25 

2.50 

3.29 

ESP* 

11.01 

2.46 

2.85 


wedge 

15.15 

2.60 

2.40 

DPP 

15.12 


asymmetric 

15.07 

3.34 

3.36 


linear* 

10.48 
rn = 2.20 
r 2 3 = 2.21 
square 
9.73 
2.33 


rhombus 

9.53 

2.40 

2.55 

pentagon 

14.25 

2.28 


linear 
13.91 
ri2 = 2.20 
r 2 3 = 2.21 

CCR 

13.09 

2.36 

3.15 

2.21 

hexagon 

18.48 

2.25 


linear 
17.34 
rn = 2.20 
r 2 3 = 2.21 
r 34 = 2.21 
PP 
16.40 
2.43 
2.53 


Td 

15.71 

2.31 

rhombus 

13.10 

2.31 

2.31 

* 

square 

10.49 

2.31 


ETB 

20.40 

2.37 

2.35 

3.84 

FSP 

20.13 

2.34 

2.32 

FTB 

16.98 

3.50 

2.34 

2.34 

octahedron 

26.52 

2.37 

FCTB* 

26.22 


OB 

23.95 


* 

square 

8.64 

2.38 

chain 
8.00 
r = 2.30 
0 = 126.75 
linear* 
7.41 
2.32 
2.34 

pentagon 

12.44 

2.32 


CCR 

11.27 

2.38 

3.50 

2.30 

FTB 

10.95 

3.36 

2.45 

3.00 

hexagon 

15.79 

2.30 

asymmetric 

14.83 

3.19 

3.26 

wedge 

14.33 

2.54 

2.43 
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TABLE IV. Equilibrium properties of optimized bulk silicon structures. AE E c Ecfdiamon ) 
and Ec is the cohesive energy (eV/atom). a (A) and B (Mbar) axe the lattice parameter and butt 
modulus, respectively, sc is die internal parameter of die BC-8 structure; it is given in units of a. In 
the DFT column, the results for BC-8 and graphite are from Refs. 101 and 102, respective y, e 
values of B for simple hexagonal and (5-tin are from Ref. 103, and the remaining data are from 
Ref. 89. The experimental values of E c and a for diamond are -4.63 eV and 5.429 A. respective y. 



•4.6297 

5.432 


Hexagonal AE 

Diamond a 
c/a 

BC-8 AE 

a 

x 

B 

[3-tin AE 

a 
c/a 
B 

Simple AE 

Hexagonal a 
c/a 
B 

Simple AE 

Cubic a 

BCC AE 

a 

HCP AE 

a 
c/a 

FCC AE 


AE 

a 

c/a 

B 


0.016 

3.858 

1.633 

0.13 

6.67 

0.1003 

0.96 

0.266 

4.822 

0.552 

1.19 

0.293 

2.639 

0.94 

1.06 

0.348 

2.528 

0.525 

3.088 

0.552 

2.735 

1.633 

0.566 

3.885 

0.710 

3.895 

2.726 

0.50 


0.0012 

3.846 

1.630 

0.311 

6.682 

0.1034 

2.73 

0.576 

5.243 

0.464 

2.71 

-0.559 

2.405 

2.866 


0.447 

2.549 

0.916 

3.165 

0.110 

4.094 

0.591 

0.950 

3.984 

0.3177 

4.123 

1.215 

2.738 


0.012 

3.841 

1.639 

0.238 

6.730 

0.1015 

2.05 

0.218 

5.113 

0.522 

3.08 

0.191 

2.762 

0.956 

3.11 

0.158 

2.609 

.312 

3.236 

0.052 

3.973 

0.685 

0.255 

4.075 

0.3847 

4.145 

1.232 

2.076 


0.0 

3.840 

1.633 

0.201 

6.591 

0.1016 

0.85 

0.213 

4.969 

0.561 

4.43 

0.403 

2.833 

0.918 

3.61 

0.293 

2.612 

0.300 

3.245 

0.321 

3.647 

0.884 

0.423 

4.147 

0.6715 

4.073 

1.193 

1.667 


0.0 

3.841 

1.633 

0.207 

6.588 

0.1054 

0.96 

0.350 

4.978 

0.521 

2.97 

0.371 

2.634 
0.997 

1.29 

0.388 

2.529 

0.479 

3.153 

0.637 

2.800 

1.633 

0.628 

3.960 

0.3625 

4.094 

1.193 

0.944 


0.0 

3.841 

1.633 

0.026 

6.579 

0.1018 

1.10 

0.455 

4.987 

0.518 

3.40 

0.527 

2.613 

0.985 

1.40 

0.343 

2.501 

0.644 

3.126 

0.551 

2.730 

1.633 

0.548 

3.861 

0.5087 

4.101 

1.198 

1.003 


0.0 

3.841 

1.633 

0.245 

6.644 

0.1008 

1.03 

0.327 

4.905 

0.524 

1.38 

0.469 

2.699 

0.967 

1.38 

0.318 

2.544 

0.432 

3.084 

0.761 

2.756 

1.633 

0.761 

3.897 

0.5131 

4.096 

1.233 

0.984 


Graphitic 

Silicon 



TABLE V. Elastic and vibrational properties of silicon. The bulk modulus and 
elastic constants are in Mbar, the phonon frequencies in THz. B' is the pressure 
derivative of the bulk modulus, C° 4 4 is the theoretical value obtained for C 4 4 m 
the absence of internal strain, and C is Kleinman's internal strain parameter; their 
value in column 2 is from Ref. 89, 110, and 112, respectively. The experimental 
values for B and the elastic constants are from Ref. 109. The phonon frequencies 
were taken from Ref. Ill, 34, 10, 27, and 13 for experiment, PTHT, BH, SW and 

T2, and T3, respectively. 


Experiment 

PTHT 

BH 

SW 

DOD 

T2 

T3 

B 

0.99 

2.788 

1.692 

1.083 

0.884 

0.98 

0.98 

B' 

4.2 

7.82 

5.66 

2.93 

4.27 

4.58 

4.30 

Ci 1 

1.67 

2.969 

2.042 

1.616 

1.206 

1.217 

1.425 

C 12 

0.65 

2.697 

1.517 

0.816 

0.722 

0.858 

0.754 

C 44 

0.81 

0.446 

0.451 

0.603 

0.659 

0.103 

0.690 

/^O 

C 44 

1.11 

2.190 

1.049 

1.172 

3.475 

0.923 

1.188 

C 

0.74 

1.03 

0.74 

0.63 

1.06 

0.83 

0.67 

VTA (X) 

4.4 

4.5 

5.6 

6.7 


2.7 

9 

VJO (X) 

13.9 

19.3 

14.5 

15.9 


15.3 

16 

VLOA (X) 

12.3 

13.8 

12.2 

13.1 


11.7 

12 

vlto CO 

15.3 

18.3 

16 

18.1 


16.5 

16 
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TABLE VI. Formation energies ( in eV) of intrinsic defects in silicon. The first 
and second values are for the equilibrium (relaxed) and ideal (unrelaxed) 
configurations, respectively. The third value is the radial relaxation of nearest- 
neighbors around the defect (in %); a negative value indicates an inward 
relaxation towards the defect. 1 T , Ih, 1b, and Is are the tetrahedral, hexagonal, 
bond-centered, and split interstitials, respectively. The DFT energies and 
relaxation for the split vacancy and the unrelaxed vacancy formation energy are 
from Ref. 96; the others are from Ref. 94. 



DFT 

PTHT 

BH 

SW 

DOD 

T2 

T3 

Vacancy 

3-4 

4.5 

0.77 

2.50 

2.12 

3.83 

2.82 

4.63 

2.57 

3.23 

2.81 

2.83 

3.70 

4.10 



38.5 

-25.7 

-24 

14.7 

1 

10.5 

Split Vac. 

4.19 

2.83 

2.30 

3.36 

4.17 

1.40 

3.50 

5.01 

4.53 

4.72 

6.00 

8.12 

4.15 

10.5 


-9.5 

-15.9 

-12.5 

-11.8 

-14.5 

-14.9 

-8.8 

It 

5 - 6 

0.63 

1.56 

5.25 

3.03 

5.03 

3.45 

a l 


1.91 

4.57 

12.21 

5.00 

5.85 

6.92 



3.8 

8 

9 

9.1 

7.3 

10.5 

Ih 

4 - 5 

0.84 

5.32 

2.89 

9.31 

6.95 

17.10 

2.61 

5.11 

3.67 

5.39 

4.61 

8.22 



7.4 

11.5 

14.7 

7.3 

7.6 

10.2 

Ib 

4 - 5 

1.92 

2.54 

5.99 

4.39 

2.84 

5.86 

is 


1.47 

3.30 

5.62 

3.49 

2.32 

4.70 
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TABLE VII. Properties of Si(100) surfaces, y is the surface energy , Ay is the relative energy 
with respect to the ideal lxl surface, and o is the lateral surface stress tensor. The x and y 
directions run along the dimer bond and the rows of dimers, respectively. Energies and stresses are 
given in eV/lxl cell. A is the first interlayer contraction (in %). rd and rbb (m A) are the bond 
lengths of the dimer and the back bond between surface and second-layer atoms. 01, 02. and 03 , 
are the bond angles as indicated in Fig. 9. The DFT results for the surface energy and the structural 
parameters of the (2x1) structure are from Ref. 90 and 125, respectively, those for the surface 
stress are all from Ref. 86, and the remaining data are from Ref. 123 (see text). 



DFT 

PTHT 

BH 

sw 

DOD 

T2 

T3 




ideal lxl 




Y 

<?xx 

Oyy 

2.5 

2.535 

0.855 

1.805 

1.176 

2.363 

2.080 

1.421 

1.683 

2.315 

0 

0 

1.779 

0 

0.145 

2.015 

0 

0.625 

2.126 

0 

-0.236 




relaxed lxl 




Ay 

^XX 

Oyy 

A 

-0.03 

-5.1 

-0.077 
-0.427 
-2.176 
- 7.0 

-0.027 

0.848 

0.273 

-5.5 

0 

0 

0 

0 

-0.085 

0.515 

-2.775 

-10.2 

-0.004 

0.023 

0.080 

-2.3 

-0.037 
0.076 
-1.693 
- 7.2 





2x1 




ay 

Oxx 

Oyy 

A 

r d 

rbb 

01 

02 

03 

-0.93 

0.693 

-1.945 

-24.4 

2.23 

2.29 

107.8 
92.9 

100.8 

-0.690 

-0.808 

-1.731 

-23.3 

2.339 

2.313 

109.0 

91.0 

104.7 

-0.709 

0.669 

0.008 

-13.3 

2.403 

2.352 

106.7 

94.8 

103.5 

-0.899 

1.167 

-0.051 

-8.3 

2.404 

2.367 

104.8 

97.9 

101.0 

-0.714 

-0.094 

-1.709 

-22.9 

2.318 

2.314 

109.0 

90.6 

106.5 

-1.258 

0.703 

0.190 

-14.6 

2.328 

2.340 

106.6 

94.8 

103.2 

-0.759 

0.367 

-1.236 

-15.6 

2.365 

2.336 

106.7 
94.7 

102.8 





c-2x2 




Ay 

Oxx 

Oyy 


-0.839 

-1.356 

-1.419 

-0.703 

0.898 

0.851 

-0.824 

1.691 

0.574 

-0.720 

0.274 

-0.866 

-1.143 

1.517 

0.567 

-0.753 

0.865 

-0.344 



Pandey 7 t-bonded defect 

structure 



Ay 

Gxx 

Oyy 

-0.895 a 


-0.687 

0.130 

2.577 

-0.814 

0.782 

3.454 

-0.045 

-1.110 

-0.227 

-1.289 

1.075 

1.441 

-0.682 

-0.061 

3.075 


a Ref. 127 


7 0 



TABLE Vin. Properties of Sid 1 1) surfaces. For the W3W3) surface. V stands for the vacancy 
model, y is the surface energy . Ay the relative energy with respect to the ideal lxl surface, and o is 
the lateral surface stress tensor. For the 2x1 it-bonded surface, the x and y directions are p e an 
perpendicular to the surface chains, respectively. Energies and stresses are given in eV/lxl ce . is 
the fust interlayer conuaction (in %). Unless indicated otherwise, the DFT results are from Ref. 8 
and 91 where the stresses were calculated at 8 Ry and the energies at 12 Ry. 


DFT PTHT 





ideal 

(lxl) 

7 

1.56 a 

0.831 

1.035 

1.158 

t 

a 


2.323 

0.969 

0 




relaxed 

(lxl) 

Ay 

-0.17 a 

-0.148 

-0.012 

0 

a 

-0.54 

-2.148 

0.149 

0 

A 

-27.0 

- 27.6 

-8.4 

0 




(2x1) 7t' 

•bonded 

Ay 

-0.45 

-0.150 


0.373 

Gy y 

1.4 

0.252 


0.017 

v AA 

Gyy 

0.4 

0.956 


2.236 




(7x7) 

DAS 

Ay < 

-0.45 

0.240 

0.398 

0.532 

g 


2.172 

1.297 

1.972 




(7x7) DS 

Ay 


-0.120 

0.093 

0.062 

G 


0.295 

0.534 

0.980 




(2x2) 

H 

> 

-<i 

-0.44 

0.085 

0.333 

0.513 

G 

1.66 

0.949 

0.606 

1.113 

H 3 Ay 

-0.33 

-0.184 

0.191 

0.267 

G 

1.18 

-3.606 

-0.397 

0.234 




(V3xV3) 

T 4 Ay 

-0.38 

0.206 

0.442 

0.456 

G 

1.70 

1.449 

0.742 

2.111 

H 3 Ay 

-0.07 b 

-0.178 

0.260 

-0.043 

G 


-4.143 

-0.502 

0.724 

V Ay 

O 
*— » 

O 

0.44 

0.45 

0.08 

G 


-2.35 

5.26 

5.98 


DOD 


0.806 

0.140 


-0.134 

-1.664 

-31.7 


0.002 

-0.770 

1.087 


0.390 

1.153 


-0.024 

-0.272 


0.320 

-0.074 

0.337 

-3.309 


0.447 

0.371 

0.478 

-3.528 

0.43 

-1.32 


0.707 

0.625 


-0.005 

0.102 

-6.4 


0.091 

0.473 

1.043 


-0.170 

1.614 


-0.052 

0.559 


-0.081 

1.248 

-0.115 

1.391 


-0.109 

1.663 

-0.143 

1.705 

0.45 

- 0.01 


1.026 

-0.074 


-0.069 

-1.241 

-20.3 


0.145 

-0.663 

1.711 


0.625 

0.738 


0.041 

0.112 


0.566 

-0.608 

0.350 

-0.621 


0.774 

-0.388 

0.482 

-0.522 

0.55 

-1.17 


a Ref. 149 

b Ref. 93, calculated at 6 Ry. 
c Ref. 148, calculated at 10.5 Ry. 



TABLE IX. Properties of the Si(110) surface, y is the surface energy , Ay the 
relative energy with respect to the ideal lxl surface, and a is the lateral su ace 
stress tensor. The x and y directions are parallel and perpendicular to t e 
surface chains, respectively. Energies and stresses are given in eV/lxl cell. A is 
the first interlayer contraction (in %). 



PTHT 

BH 

SW 

DOD 

T2 

T3 




ideal 

(lxl) 



7 

0.885 

1.043 

1.158 

0.806 

0.707 

1.026 

i 

Oyy 

2.170 

0.934 

0 

0.140 

0.625 

-0.074 

v AA 

Oyy 

1.273 

0.714 

0 

0.385 

0.468 

0.179 




relaxed 

(lxl) 



Ay 

-0.080 

-0.018 

0 

-0.068 

-0.009 

-0.035 

Oyy 

-1.627 

0.063 

0 

-1.360 

0.018 

-1.027 

v AA 

Ovrv 

-1.091 

-0.066 

0 

-0.634 

0.063 

-0.612 

A 

-5.3 

-3.4 

0 

-6.6 

-2.4 

-4.3 
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Figure Captions 


Figure L 
Figure 2. 
Figure 3. 

Figure 4. 
Figure 5. 


Figure 6. 
Figure 7. 
Figure 8. 
Figure 9. 


Illustration of the geometry of a triplet of atoms used in the definition 
of the three-body potentials. 

Comparison of the two-body potential functions, V 2 (r). The open 
circles correspond to the ab initio calculation of Ref. 83. 

Comparison of the angular variation of the three-body potentials, g(9), 
(an isoceles triangle was used for PTHT). 

Three-body energy (in eV/atom) versus angle 6 for a triplet of atoms 
forming an isoceles triangle. 0 is the apex bond angle between the 
two equal bond lengths fixed at the equilibrium bond length of the 
diamond structure of silicon, (see text) 

Illustration of the geometry of some structures for Si 5 - Si 6 . These 
are: the rhombus (1), the corner-capped triangle (2), the chain (3), 
the D2d structure of Ref. 88 (4), the flat trigonal bipyramid (5), the 
comer- (6) and edge- (7) capped rhombus, the square pyramid (8), 
the edge- (9) and face- (10) capped trigonal bipyramid, an 
asymmetric structure (11), a distorted (12) and regular (13) 
pentagonal pyramid, the orthorhombic bipyramid (14), and the wedge 
or trigonal prism (15). In the elongated form of the trigonal bipyramid, 
atoms 1, 2, and 3 are bonded to each other and the apex atoms are 
widely spaced. 

Binding energy per atom versus the number of atoms in the cluster in 
Si 2 - Si 6 . The curve KR corresponds to the ab initio results of Ref. 88. 

Comparison of cohesive energy for various bulk silicon structures. 

For the graphitic and HCP structures, see text. 

Comparison of energy versus first-nearest neighbor distance for the 
cubic diamond structure. 

Side view of the Si(100)-2xl surface showing the symmetric dimer 
reconstruction. The arrows indicate the atomic displacements (in A) 
from the bulk terminated positions. The four values (starting from the 
top) correspond to DFT(Ref. 125), BH, SW, and T3, respectively. The 
displacements for T2 are similar to those for T3; those for PTHT and 
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DOD are somewhat different. The cartesian coordinate system is 
also shown; the y direction runs into the page and is along [TlO]. 

Figure 10. Side view of the Pandey 7 i-bonded defect structure for the Si(100) 
surface. The cartesian coordinate system is the same as in Fig. 9. 

Figure 11. Relative energy with respect to the ideal lxl surface and surface 
stress for the Si(l 11) (2n+l)x(2n+l) DAS (open circles) and 
DS(filled circles) structures as a function of n for the T2 potential. 
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Figure 6 
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Figure 7 


81 




Figure 8 
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Figure 10 
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